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A BSTBACT 


This  thesis  treats  the  problem  of  determining  the  eleva¬ 
tion  angle  needed  to  initialize  an  underwater  sound  ray 
tracing  algorithm  used  to  locate  tne  position  af  a  target 
vehicle.  At  regularly  spaced  time  intervals  the  vehicle 
pings  a  synchronized  sound  signal  waich  is  received  by  a 
(short  base  line)  sonar  array  containing  four  hydrophones 
positioned  at  four  of  the  corners  of  a  cube.  Tie  wavefront 
direction  angles  are  determined  fro m  the  arrival  times  at 
tne  four  hydrophones. 

Current  methods  for  using  such  time  data  t?  produce  an 
apparent  position  suitable  for  ray  tracing  are  reviewed. 
Than  four  new  methods  are  developed  and  documented  mathemat¬ 
ically.  All  methods  are  compared  under  a  simulated  environ¬ 
ment  of  a  sound  speed  profile  which  is  linear  with  depth. 
One  of  the  new  methods  is  judged  to  be  an  impri  vemer.t  over 
current  methods  in  this  idealized  environment.  Finally  the 
improved  method  is  used  to  estimate  the  variability  in  the 
time  data  from  a  real  hydrophonic  tracking  problem. 
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I.  INTRODUCTION  AND  BACKGROUND 

A-  BACKGROUND 

Suppose  that  an  underwater  acoustical  source  emits  a 
signal  at  a  known  time.  Suppose  that  the  times  of  arrival 
of  that  signal  at  the  hydrophones  on  a  three  dimensional 
sensing  array  are  also  known.  Finally  suppose  that  the 
speed  of  sound  in  water  is  modelled  as  being  homogeneous 
over  time  and  horizontal  displacement,  varying  only  with 
depth.  Then  the  angle  of  elevation  (A),  and  the  time  (T)  of 
the  arrival  of  the  signal  at  the  acoustic  center  of  the 
hydrophone  array  can  be  determined.  If  the  relationship 
between  the  speed  of  sound  and  the  depth  under  water  is 
known  exactly,  then  the  angle  A  and  the  time  T  can  be  used 
to  trace  the  signal  trajectory  back  over  its  ray  patn 
(called  ray  tracing)  to  determine  the  original  position  of 
the  acoustical  source  [Ref.  1].  However,  full  realization 
of  this  method  in  actual  hydrophonic  tracking  ls  prevented 
by  two  primary  sources  ox  inaccuracies  in  the  process.  The 
first  and  probably  greatest  problem  is  that  tie  speed  of 
sound  profile  can  be  approximated  at  only  a  few  locations, 
usually  at  great  cost  to  achieve  even  moderate  accuracy.  In 
addition  the  profile  is  certain  to  flucuate  over  time  and 
location.  The  second  problem,  confounding  the  first,  is 
that  there  may  be  innacurac ies ,  of  unknown  size,  in  the  time 
data  values  recorded  by  the  sensing  array. 

B.  PURPOSE 

As  noted  in  [Ref.  1]  the  procedure  of  de:ermiaing  a 
sound  source  position  by  ray  tracing  is  very  sensitive  to 
even  small  errors  in  the  an  gle  of  elevation  or  speed  of 
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sound  at  the  sensing  array.  Cf  course  ray  trar  ing  is  also 
highly  iependent  on  the  accurate  de  ter  ir.inat  ion  or  the  the 
transit  time  from  source  to  array.  The  ^urpase  of  this 
stud}  is  to  develop  an  appropriate  model  of  the  timing  data 
observed  in  the  hydrophonic  tracking  problem.  Tie  objective 
of  the  desired  model  is  to  estimate  the  error  in  the  time 
values,  and  produce  improved  estimates  of  the  initial  angle 
and  ray  transit  time,  so  as  to  reduce  tae  effect  of  these 
errors  on  the  ray  tracing  procedure. 

In  pursuit  of  these  objectives  this  thesis  first  reviews 
the  currently  used  models,  which  are  called  the  NAVY  and 
NAY Y_ A  methods.  TheE  four  alternative  models  are  developed, 
called  the  L.S.,  L.S.C.,  M.  1.  ?.  and  M.L.S.  methols.  Finally 
the  performance  of  all  six  models  are  evaluated  and  compared 
using  simulation  studies.  The  comparisons  are  made  under 
the  idealized  conditions  of  a  known  linear  spea d  of  sound 
versus  depth  relationship. 

C.  TEACHING  RANGE  CCNFIG3E  ATI ON 

The  tracking  range  which  supplied  data  for  this  study 
consists  of  several  separate  three  dimensional  hydropnone 
arrays  sitting  on  the  sea  bottom.  They  are  dared  out  in  a 
rough  grid,  each  array  being  approximately  750)  feet  from 
the  next,  so  that  the  sound  source  being  tracked  is  never 
more  than  about  5000  feet  from  the  nearest  array.  The 
arrays  are  at  depths  cf  roughly  1000  to  1300  feet. 

Each  array  (see  figure  1.1)  has  four  independent  hydro¬ 
phones  defining  an  orthogonal  coordinate  system.  The  phones 
are  referred  to  as  the  x,y,  z,  and  c  hydrophones,  and  are  on 
four  adjacent  vertices  of  a  cube  (see  figure  1.2  :  with  sides 
of  length  D  (usually  30  feet).  The  arrays  are  linked  to 
shore  based  computers  by  electronic  cade.  Tm  origin  of 
the  local  coordinate  system  of  each  array  is  at  the  acoustic 
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Figure  1.1  Acoustic  Signal  Detection  by 
Three  Dimensional  Hydrophonic  Array. 
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Figure  1.2  Geometry  of  Hydrophonic  Arnys. 


The  sound  source  is  equipped  with  a  clock  synchronized 
with  the  shore  based  computers  and  emits  a  signal  at  speci¬ 
fied  intervals.  The  times  of  receipt  of  those  signals  at 
the  four  hydrophones  are  recorded,  and  the  corresponding 
travel  times  are  calculated  by  subtraction  of  the  signal 
emission  time. 

D.  APPAEENT  POSITIONS 


The  first  step  in  estimating  the  position  of  a  sound 
source  is  to  do  so  under  the  assumption  that  thi  sound  wave 
travelled  its  entire  trajectory  through  water  which  had  a 
copctdiit  speed  of  sound.  The  result,  called  the  apparent 
position,  is  obviously  erroneous,  but  is  then  corrected  by 


the  ray  tracing  procedure  (a  description  of  w: ich  follows 
later).  The  constant  velocity  value  used  is  usually  the 
estimated  speed  of  sound  V  at  the  depth  of  ;  he  sensing 
array. 

The  constant  velocity  assumption  implies  that  the  wave- 
front  is  an  expanding  sphere.  Therefore  the  apparent  posi¬ 
tion  is  calculated  using  simple  spherical  equations 
involving  squared  distance  calculations.  Specifically,  if 
Tx  is  the  travel  time  recorded  for  the  x-phone,  then  7«Tx 
is  the  distance  between  the  apparent  position  and  the 
x-phone.  This  distance  is  also  equal  to  the  usual  geometric 
distance  between  the  two  positions 

(  X  ,  Y  ,  Z  )  (apparent  position) 

and  (  D  ,-D  ,-D  )  /  2  (x  hydrophone  position). 
Equating  these  two  squared  distances,  and  equating  their 
counterparts  for  the  other  three  hydrophones,  tie  equations 
(1.1)  are  obtained. 


(1.1) 


(  X  -  D/2  )  +  (  Y  +  D/2  )2  +  (  Z  +  D/2  ) 

(  X  +  D/2  )2  +  (  Y  -  D/2  )2  +  (  Z  +  D/2  )' 

(  X  +  D/2  )2  +  (  Y  +  D/2  )2  +  (  Z  -  d/2  )' 

(  X  +  D/2  )2  +  (  Y  +  D/2  )2  +  (  Z  +  D/2  )' 


Assuming  that  the  times  Tx,  Ty,  Tz,  and  I;  are  known, 
(1.1)  is  a  system  of  four  equations  in  three  unknowns  X,  Y, 
and  Z.  This  overdetermined  system  will,  in  general,  have  no 
exact  solution.  In  fact,  even  if  the  time  values  were 
exactly  correct,  the  system  would  still  have  no  exact  solu¬ 
tion.  This  is  because  the  equations  correspond  to  the 
straight  line  ray  paths  due  to  the  constant  velocity  assump¬ 
tion,  whereas  the  time  values  correspond  to  tie  true  ray 
paths  which  are  not  straight  due  to  the  actual  variation  of 
velocity  of  sound  along  the  ray  path.  This  is  a  subtle,  but 
very  important  point. 


so  as  to 


l 


i 


To  throw  out  one  of  tne  equations  arbitrarily, 
reduce  the  system  to  three  equations  in  three  unknowns,  is 
to  throw  away  information  from  one  of  the  hydrophones.  The 
psuedo  solution  currently  utilized  is  to  subtract  the  fourth 
equation  from  each  of  the  first  three,  yielding  a  system  of 
tnree  equations  in  three  unknowns  waich  allots  an  exact 
solution  involving  information  from  all  four  hydrophones. 
However  that  solution  will  not,  in  general,  satisfy  any  of 
the  original  four  equations,  and  is  only  one  of  many  reaon- 
alle  ways  to  choose  an  approximate  solution. 


E.  INITIAL  ELEVATION  ANGLE  AND  RAY  TRANSIT  TIME 


Assuming  that  a  solution  (Xa,Ya,Za)  has  beei  determined 
for  the  apparent  position,  then  the  initial  angle  of  eleva¬ 
tion  is  just  the  angle  of  elevation  of  that  solution,  given 
by  (1.2). 


=  arcsin 


(1.2) 


The  objective  is  to  find  an  apparent  positioi  (Xa,Ya,Za) 
such  that  (1.2)  computes  an  angle  which  appraximates  the 
physically  correct  elevation  angle  as  closely  as  possible. 
The  solution  (Xa,Ya,Za)  and  the  times  lx,  Ty,  Tz  and  Tc  car. 
then  be  used  to  determine  an  appropriate  value  f> r  T,  whicn 
is  the  'ray  transit  time',  or  time  of  arrival  of  the  sound 
wave  at  the  acoustic  center  of  the  sensing  array.  The 
currently  employed  method  uses  the  proportional  relationship 
of  equation  (1.3),  where  S  and  Rc  are  the  ran)  es  from  the 
apparent  position  to  the  acoustic  center  and  to  the  c  hydro¬ 
phone  respectively,  as  in  (1.4). 
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-yf(Xa+ D/2)  +  (Va+D/2)  + 


(2+0/2) 

a 


F.  RAY  TRACIBG- 

Whichever  method  is  selected  to  produce  tie  apparent 
position  (Xa,Ya,Za),  it  is  transformed  into  the  estimate  of 
the  true  sound  source  position  (X,  Y,  Z)  by  the  procedure  of 
ray  tracing.  When  there  is  velocity  layering  ii  the  water, 
the  ray  path  is  no  longer  a  straight  line.  This  is  treated 
using  repeated  applications  of  Snell’s  Law  [Ref.  2  p.131], 
starting  with  the  layer  of  water  in  which  the  array  sits, 
and  backtracking  upwards  through  successive  velo:  ity  layers, 
until  tne  estimated  ray  transit  time  T  is  consumed. 

The  layering  effect  is  artificially  induced  by  the  limi¬ 
tation  that  the  speed  of  sound  can  be  estimated  at  only  a 
finite  number  of  depths,  the  result  of  which  is  commonly 
called  the  water  column.  For  example,  at  the  tracking  range 
studied  the  speed  of  sound  is  measured  every  25  feet, 
starting  at  tne  depth  of  12.5  feet.  Hence,  for  example,  the 
third  layer  from  50  to  75  feet  deep  is  assumed  to  have  a 
constant  speed  of  sound  equal  to  that  measured  at  62.5  feet. 

The  first  layer  processed  [Ref.  3  p. 4 ]  is  the  partial 
layer  lying  between  the  array  and  the  deepest  laf  er  boundary 
that  is  shallower  than  the  array,  with  thickness  Z1  (see 
figure  1.3). 

The  incremental  slant  range  in  the  first  Layer  is  31 
given  by  (1.5),  where  A1  is  the  initial  elevation  angle 
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Figure  1.3  Bay  Tracing  an  Apparent  Position. 


Sl  =  hi  sin  (V  (1.5) 

The  incremental  travel  time  in  tne  first  Layer  is  T1 
given  by  (1.6),  where  VI  is  the  velocity  estimated  for  the 
layer  in  which  the  array  is  situated. 


* 


in 


The  incremental  horizontal  distance  travellei 
the  first  layer  is  Hi  3iven  by  (1.7). 


(I.o) 

by  the  ray 


H1  =  s1  cos  (A_,  ) 


(1.7) 


To  determine  the  angle  of  elevation  in  the  next  layer, 
Snell’s  Law  (1.8)  is  applied,  where  7  2  is  the  speed  cf  sound 
estimated  for  that  layer. 


cost.^)  cos(A,) 


When  (1.8)  is  solved 
entry  into  the  next  layer, 

cos (A,) 


for  the  cosine  of  tie  angle  of 
(1.9)  is  obtained. 
v0  cos  (A.,  ) 


The  procedure  of  computing  the  incremental  values  of 
slant  range,  time  and  horizontal  distance  are  now  repeated 
for  the  second  layer.  The  overall  procedure  is  repeated 
upwards  through  layers  2,...,n  ,  where  n  is  the  first  layer 
in  which  the  sum  of  the  incremental  travel  times  exceeds  tne 
total  ray  transit  time,  as  in  (1.10). 


+  .  .  .  +  >  T 


(1. 10) 


In  the  last,  uppermost  layer  the  values  Tn,  Sn,  Hn,  and 
Zn  must  he  adjusted  to  compensate  for  overshooting  the  total 
time  T.  The  values  Hi  and  Zi  (i=1,...,n)  are  then  accumu¬ 
lated  to  get  (1 .  11)  . 
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(i.  ii) 


Now  the  ray  traced  position  estimate  is  given  by  (1.12). 
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not  due  to  the  approximation  by  layers,  except  possibly  when 
there  are  radical  changes  in  the  velocity  pattern  within 
single  layers  such  as  frequently  occur  in  layers  near  the 
water  surface.  Rather  such  errors  apparently  ara  due  mostly 
to  inaccuracies  either  in  the  estimation  of  the  speed  of 
sound  profile  itself,  or  in  the  initial  angle  and  transit 
time  estimates.  This  study  will  focus  on  those  arrors  which 
are  involved  in  the  time  values  observed  at  the  hydrophones, 
and  attempt  to  produce  methods  of  initial  angle  and  transit 
time  estimation  which  reduce  tne  effects  of  tnoae  errors  on 
the  overall  position  estimation  procedure. 
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II.  COBB  ENT  LI  EMPLOYED  METHODS 


A.  EASIC  CONCEPTS 

The  met hoi  currently  used  for  estimation  of  an  initial 
angle  and  ray  path  transit  time  focuses  or.  the  four  spher¬ 
ical  equations  (1.1)  given  in  Chapter  I.  A;  previously 
discussed,  these  equations  have  no  exact  solution  because: 

1.  they  form  an  overdet eroined  system  of  four  equations 
in  three  unknowns; 

2.  the  time  values  recorded  at  the  hydrophones  may  be 
innacurate  due  to  unknown  sources  of  error  in  the 
observation  process;  and 

3.  even  if  the  time  values  were  exact,  they  correspond 
to  a  nonconstant  velocity  profile,  and  so  will  not  be 
be  correct  for  the  constant  velocity  geometry  (spher¬ 
ical  wavefront)  assumed  by  the  equations. 

As  previously  mentioned,  the  psuedo  solatia  n  chosen  by 
the  current  method  is  to  subtract  the  fourth  spherical  equa¬ 
tion  from  each  of  the  first  three,  and  solve  the  resulting 
system  of  three  equations  in  tnree  unknowns.  This  method 
has  the  beneficial  quality  that  information  is  retained  from 
all  four  hydrophones,  whereas  to  just  drop  the  fourth  equa¬ 
tion  (or  any  one  of  the  equations)  without  the  initial 
subtraction  would  cause  complete  loss  of  the  information 
from  the  data  recorded  by  one  of  the  hydrophones.  However 
it  is  important  to  realize  that  the  solution  thus  obtained 
does  not  actually  satisfy  any  of  the  original  four 
eguat ions. 

It  should  be  noted  that  the  development  of  ti  is  solution 
in  [fief.  3]  is  done  entirely  from  a  geometrical  point  of 
view,  and  does  not  mention  the  system  of  foi  r  spherical 


equations.  The  text  of  [Ref.  3]  Joes  not  draw  i  ttentior.  to 
the  fact  that  the  solution  developed  is  just  one  of  many 
plausible  choices,  none  of  which  will  satisfy  all  four 
spherical  constraints.  Therefore  the  solutioi  chosen  is 
treated  as  though  it  were  the  exact  solution,  only  subject 
to  errors  in  the  observed  hydrophone  time  values.  However, 
even  with  exactly  correct  time  values,  tnis  currently 
employed  method  will  not  yield  the  true  elevation  angle  and 
ray  transit  time.  This  is  due  to  the  assumption  of  a 
constant  velocity  vice  the  true  nonconstant  velocity 
profile.  This  conflict  introduces  an  automatic  bias  in  the 
initial  angle  and  time  estimates  currently  used  for  ray 
tracing. 


B.  CCHPOTATIONS 


To  simplify  notation,  let  (X,Y,Z)  be  the  coordinates  of 
the  apparent  position  which  was  formerly  denoted  (Xa,Ya,Za). 
Then  when  the  current  solution  is  applied,  the  first  step  is 
to  subtract  the  fourth  equation  from  eacn  of  the  other 
three,  which  produces  the  equations  (2.1). 


(  X  -  D/2  )' 


(  X  +  D/2  )' 


=  V“  ( 
2  . 


(  Y  -  D/2  )  -  (  Y  +  D/2  )’  =  V 

2  .  _  .  ..  .2 


„2 

‘c 

2 


) 


=  ->’ 
2  2  0 
=  V  (  T  -  T  ) 
c  z 


(  Z  -  D/2  )"  -  (  Z  +  D/2  ) 

The  solution  to  these  are  easily  obtained,  as  in 


x  =  v2  (  t;  -  T2  )  /  2  D 

Y  =  V2  (  T2  -  )  /  2  D 

2  =  "2  <  -  T2  )  /  2  D 


(2.1) 


(2.2). 

(2.2) 


Then  the  initial  elevation  anale  estimate  is  (2.3). 


The  ray  patn  transit  time  to  the  acousti:  center  is 
(2.4),  wnere  ?.  and  Rc  are  as  defined  in  Cnapter  I  by  (l.u). 


T  =  -c  R  ‘--c  (2.4) 

This  method  -of  determining  the  apparent  position  sr.ali 
hereafter  be  referred  to  as  the  ’Navy  hetnod',  or  'NAVY'  for 
short . 


C.  ADJUSTMENTS  TO  THE  ORIGINAL  SOLUTION 

Experience  has  shown  that  the  NAVY  method  produces  an 
apparent  position  estimate  which  usually  can  be  improved  by 
an  adjustment  which  is  described  in  this  section.. 

The  cosine  of  the  angle  between  the  i-th  axis  and  the 
straight  line  from  the  origin  out  to  the  apparent  position 
is  called  the  i-th  direction  cosine  Ci.  It  is  a  fact  of 
geometry  that  the  sum  of  the  sguares  of  the  three  direction 
cosines  must  egual  unity.  Therefore  the  method  is  now 
adjusted  to  reflect  that  constraint. 

The  direction  cosines  used  are  the  angles  nade  by  the 
ray  path  at  the  c  hydrophone.  Therefore  the  (X,Y,Z)  coordi¬ 
nates  calculated  by  the  original  method  are  first  translated 
to  coordinates  referenced  to  the  c-pnone  as  the  temporary 
origin,  as  in  (2.5). 


X  =  X  +  D/2  Y  =  Y  +  D/2  Z  =  Z  +  D/2 

c  c  c 

Therefore  the  three  direction  cosines  ar? 


(2.5) 

given  by 


(2.  6)  . 


c 


x 


X  /VT  C  =  Y  /  V  T  C 

c  c  y  c  c  z 


Z  /  V  T 
c  c 


(2.6) 


The  denominators  in  (2.6)  are  all  V*Tc  becj  use  that  is 
the  range  from  the  apparent  position  to  the  c  hydrophone,  as 
estimated  by  the  time  from  the  c  hydrophone.  Ideally  tr.ese 


cosines  sacUiU  aca  to  unity  when  sguarea.  Inera lore  ir 


is  the  ’direction  cosines  correction'  factor  defined 


0  7 


(2-7)  ,  then  DCC  should  be  close  to  one. 

ccc  =  /c2  +  e1  7  c2 

Y  X  y  Z 


(2.7) 


deviation  of  DCC  from  unity  is  interpreted  as  an  indica¬ 
tion  of  receiver  timing  errors,  array  malformation  or 
invalid  data  at  one  or  more  of  the  hyiropnoaes  [Ref.  3 
p.C-3].  Currently  if  DCC  lies  outside  tie  interval 
(0.  98,1.02),  the  data  is  discarded  as  being  excessively  fall 
of  error.  The  direction  cosines  of  the  remaining  acceptable 
data  points  are  rescaled  using  (2.8)  to  assure  satisfaction 
of  the  direction  cosines  constraint. 


c  /  DCC 


C'  =  C  /  DC 

y  y 


c*  =  c  /  DCC 

2  Z 


(2.8) 


A  corrected  set  of  new  coordinates  are  ;  omputed  by 
(2.  9)  ,  still  being  referenced  to  the  c  nydrophone. 


.  ^  t 

‘c  "x 


y  =  V  T  C 
c  c  v 


These  are  then  translated  by 
referenced  to  the  acoustic  center. 


(2.10) 


-2  (2.9) 

to  coordinates 


x  =  x  -  d/: 

c 


V  -  D/2 
c 


y: 


(2.  10) 


This  adjusted  method  of  determining  the  apparent  posi¬ 
tion  shall  hereafter  be  referred  to  as  the  'Navy  Adjusted 
Method  '  ,  or  '  N  A  V  Y_  A  '  for  short. 


D.  DIECUSSIOH 

When  a  sound  source  is  within  the  detection  range  of 
more  than  cne  sensing  array,  each  array  produces  timing 
data.  The  data  from  each  array  can  be  processed  to  produce 


a  position  estimate.  ueanv  those  muxtip'le  estimates  o; 
position  will  he  in  reasonably  close  agreement.  However 
experience  with  actual  tracking  data  has  shown  that  this  is 
not  the  case  in  many  of  the  multiple  detection  opportuni¬ 
ties.  This  tendency  toward  disagreement  between  multiple 
estimates  ox  the  same  position  is  commonly  called  the  cross¬ 
over,  or  crosstalk,  problem.  This  problem  often  occurs  when 
the  sound  source  is  moving  away  from  the  tracking  domain  of 
one  array  into  the  tracking  domain  of  another.  This  study- 
focuses  on  improvement  of  the  initial  angle  and  time  esti¬ 
mates,  wnich  hopefully  will  help  alleviate  the  crossover 
problem. 

The  current  chcice  of  a  'best'  compromise  solution 
appears  to  be  based  on  reasons  of  simplifying  geometry  and 
calculations.  These  are  worthwhile  onjectives,  nut  do  not 
in  themselves  reflect  the  need  to  estimate  accurately  the 
initial  elevation  angle  and  ray  transit  time.  Since  there 
exist  physically  correct  values  for  both  the  ai  gle  and  the 
time,  those  values  will  produce  the  exact  position  after  ray 
tracing,  provided  that  the  velocity  profile  is  known 
exactly.  The  desire  then  is  to  estimate  these  true  values 
as  accurately  as  possible. 

The  guestion  at  this  point  is  whether  or  not  the  direc¬ 
tion  cosines  adjustment  causes  the  estimated  apparent  posi¬ 
tion  to  be  closer  to  the  true  apparent  position.  Experience 
has  indicated  that  it  does  [Ref.  3  p.C-7].  However  the 
effect  of  the  adjustment  can  be  interpreted  in  terms  of  the 
original  four  spherical  equations  (1.1).  The  rescaling  of 
the  direction  cosines  given  by  (2.6),  so  as  to  assure  that 
their  sguares  add  to  unity,  is  equivalent  to  rescaling  the 
quantities  in  (2.5)  so  as  to  assure  that  their  sguares  add 
to  ( V  «Tc) 2 .  That  is  exactly  the  constraint  stated  by  the 
fourth  spherical  equation  of  (1.1)  .  So  the  effect  of  the 
adjustment  is  to  require  that  tne  fourti  equation. 


concerning  the  data  at  hydrophone  c,  is  exactly  satisfied. 
This  requirement  will,  in  general#  assure  that  the  ether 
three  equations  are  net  satisfied.  Since  experience  shows 
that  the  adjustment  often  improves  the  solution#  this  seems 
to  imply  that  the  fourth  equation  is  somehow  more  important 
than  the  other  three.  Or  it  may  just  Le  that  exact  satis¬ 
faction  of  one  of  the  equations  usually  assures  a  reasonably 
good  compromise  solution. 

Tc  summarize#  the  NAVY  metnod  provides  a  useable 
apparent  position  suitarle  as  input  for  ray  tricing.  But 
the  direction  cosines  adjustment  used  in  the  HAVY_A  method# 
for  reasons  not  understood  at  this  time#  appears  to  improve 
that  position  as  indicated  by  test  results.  Ti ose  results 
are  supported  by  the  results  of  this  thesis  (see  Chapter  V)  . 
However,  as  will  be  demonstrated  by  the  example  considered 
in  the  next  section#  the  DCC  correction  factor  o:  the  NAVY_A 
method  has  an  effect  which  must  be  something  more  than  just 
the  smoothing  of  timing  errors. 


E.  A  COMPUTATIONAL  EXAMPLE 
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Figure  *..1  shows  the  estimated  sound  velotity  profile 
which  was  estimated  for  the  data  used  in  the  course  of  this 
thesis.  As  can  be  seen,  the  proriie  is  primarily  linear  at 
depths  greater  than  ICO  feet.  The  profile  at  deptns  greater 
than  100  feet  is  reasonably  approximated  by  the  linear 
relat ionship 

V  =  4840.7  ♦  0.03314  •  DEPTH  . 


Therefore  suppose  that  in  tnis  example  problem  the 
velocity  profile  is  known  exactly,  and  is  given  ay  the  amove 
linear  relationship. 

Under  these  circumstances,  witn  known  linear  velocity 
profile,  and  known  sound  source  location,  the  exact  times  of 
arrival  of  the  sound  wave  at  the  four  hydrophones  can  be 
computed  using  the  methods  set  forth  in  Appendix  A.  Those 
exact  times  (in  seconds)  are: 

Tx  :  0.  6779686893  Ty  :  0.674  2257788 
Tz  :  0.6782243197  Tc  :  0.6798324156  . 

The  corresponding  exact  values  for  the  initial  elevation 
angle,  ray  transit  time  and  resulting  apparent  position  are 
also  directly  computable.  Those  computations  will  hereafter 
be  be  known  as  the  EXACT  method,  and  will  produce  the 
correct  true  position  after  ray  tracing.  The  results  of  the 
EXACT  method  are  given  in  Table  I,  along  with  the  corre¬ 
sponding  apparent  position  estimates  produced  waen  the  two 
methods,  NAVY  and  NAVY_A,  are  applied  to  the  (errorless) 
time  values. 

At  first  glance  the  differences  in  Table  I  mignt  seem 
rather  small.  However  it  is  important  to  recall  that  these 
are  produced  under  the  ideal  conditions  of  a  very  smooth  and 
exactly  known  velocity  profile.  These  idealizations  are  far 
from  the  realities  cf  a  nonlinear  velocity  profile  which  is 
es  mated  by  a  procedure  involving  errors  which  are  unknown 
and  probably  significant.  Such  realities  sigh:  well  cause 
the  differences  in  Table  I  to  be  si  gmf  ica  ntly  larger.  The 
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TABLE  I 

Single  Example  Comparison  of 


NAVI  and  NAVY  A  Methods 


Method 
NAVI 
N  A  V Y_A 
EXACT 


Transit  Time 
T  Isecs)  ~ 

0. 67526969 

0. 67527027 

0. 67527043 


Slev.  Angle 

I~T2egsJ — 

15. 26459 

15. 26461 

15.  27  002 


NAVI 
N  A  VI_ A 
EXACT 


Apparent  Position  Estimate 
(  1005.  957  ,  3017.874  ,  868.  143  ) 

(  1006.087  ,  30  18.259  ,  868.25  5  ) 

(  1005.  966  ,  30  17.899  ,  868.555  ) 


nature  and  size  of  those  differences  remain  difficult  to 
determine  until  more  is'known  about  the  the  velocity  profile 
estimation  errors  and  their  effect  on  the  position  esti¬ 
mating  process.  In  any  event  even  the  small  differences  in 
Table  I  might  be  magnified  during  the  ray  tracing  process 
under  the  conditions  of  a  realistic  velocity  profile. 

The  differences  illustrate  the  very  important  point  that 
the  directi  on  cosines  adjustment  causes  changes  in  the  esti¬ 


mates  even  when  the  time  data  is  free  of  all  error. 


Hence 


the  deviation  from  1.0  of  the  correction  factor  DCC  is  not 
just  due  to  array  malfunction,  receiver  timiig  error  or 
other  sources  of  non-valid  data  as  previously  assumed. 

In  the  example  above,  the  NAVI_A  method  produced  a 
slightly  better  time  and  angle  value  than  the  1  AVI  method. 
However,  in  this  same  example  the  NAVI_A  method  produced  an 
apparent  position  estimate  which  is  sligntly  farther  away 
from  the  EXACT  answer  tuan  the  position  estimated  by  the 


NAVY  method.  This  is  only  one  example,  and  its  results 
should  not  be  generalized.  However  it  illustrates  the  point 
tnat  apparently  the  true  effect  of  tne  adjustment  may  not  be 
well  understood. 

Furthermore,  the  adjustment  seems  to  place  heavy 
emphasis  on  the  time  value  recorded  at  hydrophone  c. 
Therefore  the  effect  of  the  adjustment  may  well  depend 
largely  on  the  accuracy  of  that  one  particular  data  value, 
which  is  a  relatively  unbalanced  dependence  in  the  presence 
of  data  errors. 


III.  PLANAB  SAVEFRONT  MODELS 


A.  TEE  PLANAR  WAVEFRONT  ASSUMPTION 

When  a  sound  wave  travels  in  constant  veloci:  y  water,  it 
expands  in  the  shape  of  a  sphere.  If  the  velocity  profile 
is  variable  instead,  but  is  reasonably  well  behaved  versus 
depth,  then  the  expanding  wave  is  a  smooth  distortion  of  a 
spherical  surface.  In  either  case,  if  the  wave  has 
travelled  a  long  distance  when  it  arrives  at  a  hydrophoric 
array,  then  that  small  piece  of  the  wavefront  t  hich  passes 
through  the  30  foot  cube  spanned  by  the  array  may  be  approx¬ 
imated  reasonably  by  a  flat  planar  surface.  This  approxima¬ 
tion  is  the  basis  of  the  planar  wavefront  models  developed 
in  this  chapter. 


B.  PLANE  EQUATIONS 

A  plane  in  space  is  fully  defined  by  idej tif ying  any 
point  (X0, Y0,Z0)  on  the  plane,  and  also  a  vector  (C1,C2,C3) 
of  unit  length  which  is  perpendicular  to  tnat  plane.  The 
vector  is  called  the  unit  normal  vector  for  that  plane. 
Then  any  point  (X,  1,2)  on  that  plane  must  satisfy  the  equa¬ 
tion  of  the  plane,  namely  (3.1). 


c 


l 


X  -  X0  )  +  c2  (  Y  -  Yo 


The  perpendicular  distance 
(X 1 ,  Y  1,  Z  1)  not  on  the  plane  is 


from  the  plane  t> 
the  absolute  value 


0  (3-1) 

any  point 
of  (3.2). 


C1  (  X1  -  Xo  >  +  C2  (  Y1  "  Y0  >  +  C3  '  21  - 


(3.2; 


3  1 


c. 


BASIC  SCDSL  EQUATIONS 


Consider  a  coordinate  system  whose  origin  13  at  hydro¬ 
phone  c,  and  whose  axes  are  aligned  with  the  three  array 
aras.  Let  Cl,  C2  and  C3  be  the  direction  cosines  on  the  X, 

Y  and  Z  axes  respectively  for  the  vector  from  tie  origin  to 

the  apparent  sound  source  position.  Tnen  (C1,C2,C3)  is 
itself  a  vector,  of  unit  length,  which  is  perpendicular  to 
the  planar  wavefront  emanating  from  the  soi  nd  source. 
Therefore  (C1,C2,C3)  can  be  used  as  the  normal  vector  for 
the  wavefront  plane. 

In  the  coordinate  system  referenced  to  the  c-phone  as 
the  origin,  the  acoustic  center  has  coordinates  (D,D,D)/2, 
where  D  is  the  length  of  an  array  arm.  When  the  soundwave 

plane  arrives  at  the  acoustic  center,  it  wiL  1  have  the 

equation  (3.3)  . 

C:  X  +  C2  Y  +  C3  *  -  D  (  Cj_  +  w2  +  C3  )  /  2  (3>3) 

The  x-pnone  has  coordinates  (D,0,0),  ana  :  he  distance 
between  it  and  the  soundwave  plane  at  the  acoustic  center  is 

(3.4),  which  then  simplifies  to  (3.5). 

Cl  (  D/2  -  D  )  +  C2  (  D/2  -  0  )  +  C3  (  D/2  -  0  )  (3>4) 


D  ( 


“I 


)  /  2 


(3.5) 


This  distance  is  measured  in  a  direction  perpendicular 
to  the  wavefront  plane,  and  so  is  measured  in  the  direction 
of  travel  of  the  soundwave.  Therefore  the  same  distance  is 
also  e  gual  to  (3.6). 


V  (  T  -  T  ) 
x 


(3.6) 
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In  (3.6)  V  is  the  velocity  of  sound  at  the  array,  Tx  is 
the  time  of  arrival  of  the  sound  wave  at  the  x-phone,  and  I 
is  the  time  of  arrival  of  the  sound  wave  at  the  acoustic 
center.  The  term  "distance1’  is  used  loosely  in  (3.5)  and 
(3.6)  ,  because  these  quantities  may  be  negative.  The  true 
distances  are  the  absolute  values  of  these  quantities. 
Since  the  next  step  is  to  equate  these  two  distal  ces,  it  is 
only  necessary  to  show  that  these  two  quantities  always  have 
the  same  sign.  There  are  two  cases  to  consider,  depending 
on  whether  the  first  component  of  tne  apparent  position  is 
positive  (X >0)  ,  or  negative  (X<0)  .  let  (X',0,0)  be  the 
intersection  of  the  X  axis  with  the  wavefront  plane  as  it 
passes  through  the  acoustic  center.  Then  (3.3)  car.  be  used 
to  solve  for  X',  namely 

X’  =  D  (Cl  +  C2  +  C3)  /  2  Cl. 

Now  consider  the  case  where  X>0.  Then  C1>0  i Iso,  and  if 

(3.6)  is  positive,  then 
Tx  >  T 

wnich  implies  that  the  wave  arrives  at  the  acoustic  center 
before  it  arrives  at  the  the  x  hydrophone.  Therefore,  since 
X>0,  the  plane  at  the  acoustic  center  will  intersect  the  X 
axis  at  a  point  beyond  the  x  hydrophone,  or  X’>D,  and  hence 

X’  >  C  =>  (Cl  +  C2  +  C3)  /  2  Cl  >1 

=  >  Cl  ♦  C2  ♦  C3  >2  Cl  (sinca  C1>0) 

=  >  -C1  +  C2  +  C  3  >  0 
=  >  (3.5)  is  positive  . 

A  parallel  argument  applies  for  tne  case  of  X>0,  thus 
establishing  that  (3.5)  and  (3.6)  always  have  the  same  sign. 
Equation  (3.7)  is  the  result  of  equating  these  two 
quantities. 


-C1  +  C2  +  C  +  (  2  V  T  /  D  )  -  (  2  V  T  /  D  ) 


0 


1 1  7  \ 


For  convenience/  let  K  =  2V/0,  and  then  apaly  the  sane 
logic  to  the  distances  of  the  y,  z,  and  c  hydrophones  from 
tae  wavefront  plane  as  it  passes  through  tee  acoustic 
center,  to  obtain  the  equations  of  [3-d)  . 


-  K  T  +  K  T 


+  C3  +  K  T 


+  K  T 


+  C.  -  K  T  +  K  T 


(3.2) 


The  system  (3.3)  is  four  equations  in  four  ui knowns,  and 
is  the  planar  model's  version  of  the  equations  (1.1).  The 
unknowns  in  (3.8)  are  Cl,  C2,  C3  and  T.  However  there  is 
the  additional  constraint  that  Cl,  C2,  and  C3  are  direction 
cosines,  and  therefore  the  direction  cosines  constraint 
(3.9)  is  a  fifth  equation,  creating  a  a  system  of  five  equa¬ 
tions  in  four  unknowns. 


2  _2 
+  S  *  C3 


(3.9) 


Generally  there  is  no  set  of  values  (C1,C2,:3,T)  which 
will  satisfy  ail  five  equations  at  once.  This  is  because  of 
the  realities  of  a  nonplanar  wavefront  and  the  presence  of 
timing  errors.  The  next  section  developes  a  method  that 
produces  set  of  values  for  the  unknowns  which  is  intended  to 
satisfy  those  equations  reasonably  well. 

D.  HIHI BIZ  ATI ON  OF  SOM  OF  SQUARED  ERRORS 

Let  Ei,  (i=1,2,3,4)  be  the  value  of  the  left  hand  side 
of  the  i-th  equation  of  (3.8).  Then  Ei  measures  the  amount 
of  error  in  the  i-th  equation  caused  by  tne  chosen  solution. 
It  is  impossible  to  have  Ei=0  for  all  i=1,2,3,4.  However 
some  compromise  may  be  made.  Specifically  the  compromise 
chosen  here  is  the  classic  minimization  of  (3. ID) ,  the  sum 
of  squared  errors. 


~1 


E  +  r 
1  “2 

minimization  is  to  be  done  subject  to  the  cirectior. 
cosine  constraint  (3.9).  Application  of  the  Lagrange  multi¬ 
plier  technique  [Ref.  5  p.  5  5  ]  calls  for  the  minimization  of 
(3.11)  over  all  possible  choices  of  Cl,  C2,  C3,  T  and 

lambda. 


*  x  U  > 


(3.10 


Taking  tne  partial  derivative  of  L  with  re  spect  to  7 
yiel  ds  (3.12). 

O  l  _  v  /  _  o  v  r*  \  ir  /  r*  _l  r*  4-  P  T  \ 


1  =  V  (-2K  E.) 

d  T  i=l  1 


•2K  (  E  +  E  +  E  +  E.  ) 
12  3  4 


(3.  12) 


-2K  T  +  K  (  T1  +  T2  +  T3  -  T4  )  j 


Equating  (3.12)  tc  zero  and  solving  for  T  yields  immedi¬ 
ately  the  appropriate  estimate  (3.13)  of  T,  the  ray  transit 
time . 


(  T  +  T_,  +  T3  -  T4  ) 


(3. 13) 


The  partial  derivative  of  L  with  respect  to  Cl  is 
(3.14). 

— — —  =  2(r  -  E  -  P  +  E  )  -  2  X  C 

v  “2  34'  AH 

dc,  (3.  14) 


(3.  14) 


4C^  +  2KT  +  K 


'  *  *l_*2-'t3_'*4  ^  X  H 


If  (3.13)  is  used  for  T  in  (3.14),  tnen  (3.15)  results. 


a.  3  r  t  n  e 


If  the  same  procedure  is  used 
lives  or  L  with  respect  to  each  of  Cl, 
are  equated  to  zero,  then  the  results  are 


par:  i al 
C 2  and  C3, 

(3.  16)  . 


den  va- 
ar.d  ail 


X'  ci 


2  K  ( 


i  =  r  , . 


io) 


In  order  to  solve  rot  the  Lagrange  multiplier  lamhda, 
square  both  sides  of  the  three  equations  of  (3.16),  and  add 
the  resultinq  equations  together.  Then  use  the  sum  of 
squares  constraint  (3.9)  ana  solve  for  Iambi  a  tc  yield 

(3.17). 


r  ? 


(3. 17) 


Substitute  (3.17)  in  (3.16)  and  simplify  to  obtain  the 
appropriate  estimate  cf  Ci,  namely  (3.18). 


c. 

1 


1,2,3 


(3. 18) 


The  choice  of  sign  in  (3.18)  is  determined  by  the  fact 
that  Ci  is  positive  if  and  only  if  the  sound  wavj  arrives  at 
the  i-th  phone  before  it  arrives  at  the  c-phone,  which  in 
turn  implies  that  (T4-Ti)>3. 


E.  T BE  LEAST  SQUARES  METHOD 

In  summary,  the  first  alternative  metnod  for  determining 
an  apparent  position  starts  with  estimation  >f  the  ray 
transit  time  to  the  acoustic  center  by  (3.19).  Then  the 
apparent  position  estimates  are  computed  using  (3.20). 


This  method  shall  be  hereafter  referred  to  as  the  'Least 
Squares  Method',  or  'L.3.  '  for  short.  The  apparent  advan¬ 
tages  of  the  L.S.  method  are  that: 

1.  all  four  hydrophone  times  have  ecual  weight  in  a 
simple  expression  for  the  ray  transit  time  T,  rather 
than  using  an  expression  so  heavily  dependent  on  Tc 
as  in  the  NAVY  and  N  AVY_A  methods; 

2.  the  differences  of  squared  time  values  i  hich  appear 
in  the  solutions  of  the  NAVY  methods  are  avoided  in 
the  L.S.  method,  thereby  lessening  tie  tendency 
toward  computational  roundoff  problems; 

3.  the  direction  cosines  already  add  to  unity  when 
squared,  requiring  no  arbitrary  adjustmen:  ;  and 

4.  computation  of  the  initial  angle  allows  cancellation 
of  several  terms,  resulting  in  the  simple  expression 
(3.21). 


arcsin  (  T  -  T  } 
4  3 


Y,  <  T.  -  T  .  )  ' 


(3. 21) 


F.  BIAS  IN  THE  LEAST  SQUARES  METHOD 

Unfortunately  a  potentially  serious  problem  exists  with 


the  L.S. 


met  hod. 


That  concerns  the  consegueices  of  the 


assumption  of  a  planar  wavefront. 


The  effect  is  difficult 


to  derive  explicitly  because  it  is  difficult  :o  detersive 
tie  true  shape  of  a  wavefront  cor responiin g  to  a  nor.constant 
velocity  profile.  However  if  it  is  assumed  tnat  a  spherical 
assumption  is  more  accurate  the  planar  assumption,  tner.  the 
bias  can  be  estimated  roughly,  and  then  subtracted  from  the 
original  L.S.  position  estimate.  This  is  still  a  difficult 
problem  because,  as  previously  discussed,  the  four  basic 
spherical  equations  themselves  have  no  exact  solution. 

Nevertheless,  as  a  rough  estimate  of  the  bias,  the 
following  procedure  is  used.  First  estimate  the  apparent 
position  by  the  L.S.  method.  Tnen  calculate  the  straight 

line  distances  from  that  position  tc  each  of  tha  four  array 
hydrophones.  Divide  those  distances  by  the  velocity  of 
sound  at  the  array  tc  obtain  tne  corresponding  times.  Use 
these  times  to  recalculate  tne  apparent  position  using  the 
L.S.  method  again.  The  difference  between  the  original  and 
recalcualted  L.S.  positions  roughly  measures  the  error  that 
would  be  made  by  the  L.S.  method  wnen  it  is  applied  to  tne 
time  values  which  correspond  to  a  spherically  spreading 
sound  wave  whose  source  is  in  the  vicinity  of  the  original 
apparent  position.  Therefore  this  difference  cai  be  used  as 
a  bias  vector  which  can  be  subtracted  from  the  original  L.S. 
solution.  This  bias  correction  is  the  basis  of  the  method 
set  forth  in  the  next  section. 

G.  THE  LEAST  SQUARES  CORRECTED  METHOD 

The  second  alternative  method  for  estimating  ar.  apparent 
position  is  as  follows: 

1.  calculate  the  apparent  position  PI  by  using  the  L.S. 
method; 

2.  calculate  the  distances  from  that  position  to  the 
four  hydrophones,  and  convert  them  to  times  by 
dividing  by  the  speed  of  sound  at  the  array; 
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use  the  new  tir.es  to  calculate  a  r.ew  apparent  posi¬ 
tion  ?2  ,  using  tre  L.S.  metnod; 

calculate  the  difference  vector  P2-P1  (see  figure 
3.1),  and  subtract  it  from  the  original  position  PI 
to  obtain  the  corrected  position  ? 

P  =  PI  -  (P2-P1)  =  2P1-P2 

finally  adjust  the  ray  transit  time  T  calculated  for 
the  original  position  ?1,  by  using  the  proportional 
transformation: 

T '  =  T  *  E  /  R 1 

where  E  is  the  range  to  the  new  position  >,  and  F.1  is 
the  range  to  the  original  position  ?1. 


Figure  3.1  Bias  Adjustment  for  the  L.S.  Method. 


This  method  shall  hereafter  be  referred  to  as  the  'Leas 
Squares  Corrected  Method',  or  L.3.C.  for  short.  The  proper 
ties  of  the  L.S.C.  method,  liie  tnose  of  the  NAfY_A  method 
are  not  well  understood  at  this  time.  It  is  offered  only  a 


an  alternative  which  may  combine  the  beneficial  properties 
of  the  L.S.  and  NAVY  methods,  namely  that  it  will: 

1.  estimate  a  ray  transit  time  value  e^ualLy  dependent 
on  all  iou.  data  times,  thereoy  smoothing  out  exces¬ 
sive  error  in  any  one  of  the  time  values;  and 

2.  reflect  the  spherical  wave  assumption,  believed  to 
provide  a  more  accurate  description  than  the  planar 
assumption. 

For  targets  at  a  range  of  3000  feet,  the  length  of  the 
error  vector  P2-P1  varies  from  0  to  as  much  as  10  or  11  feet 
(see  Table  II).  The  error  vector  lengtns  seem  to  be  depen¬ 
dent  cn  both  azimuth  and  elevation  angles  of  the  target  from 
the  array.  These  patterns  indicate  a  potential  for  further 
investigation  to  relate  estimation  errors  to  suci  variables. 

H.  HAXIHUH  LIKELI HOCE  CALCULATIONS 

Cne  shortcoming  of  ail  methods  described  thus  far  is 
that  none  of  them  can  be  used  to  produce  an  estimate  of  the 
underlying  error  in  the  time  data  values.  The  method  set 
forth  in  this  section  is  a  first  attempt  to  estimate  that 
noise,  and  is  again  based  on  the  assumption  of  a  planar 
wavefront. 

Ler  Ti  be  the  time  recorded  from  by  i-th  hydrophone. 
Let  'Ji  be  tne  true  time  which,  under  absolutely  error  free 
conditions,  would  have  been  recorded  at  tne  i-th  hydrophone. 
An  assumption  of  laussian  noise  is  now  made,  namely  tnat 

Ti  =  Ui  ♦  Ei 

where  the  Ei  are  independent  identically  distributed  normal 
random  variables  with  mean  zero  and  variance  32.  Therefore 
the  Ti  (i=1, 2, 3, 4)  are  also  normally  distributed  with  the 
same  variance,  but  with  means  Ui  (i  =  1,2,3,4). 


TABLE  II 

Error  Vector  Lengths  for  the  L.S.C.  Method 
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Now  let  Cj  be  the  j-th  direction  cosine  (j=1,2,3).  Then 
geometrically,  using  the  assumption  of  a  planar  wavefront, 
Cj  is  given  by  (3.22),  where  V  is  tae  speed  cf  sound  at  the 
array,  and  D  is  the  array  dimension  (30  feet)  . 


v  (  u .  -  u .  )  '  D 
-*  : 


(3.  22) 


Letting  U  =  04,  then  (3.22)  can  be  solved  for  each  Uj  in 
terms  of  U  and  Cj,  as  in  (3.23). 


U  .  =  U-DC./V 

3  3  ' 

The  probability  density  of  each  time  value 
by  (3.24). 


(3.23) 


Ti  is  given 


w 


V^T 


-  (  T.  -  U.  )' 

l  l 


(3.24) 


Using  (3.23)  in  (3.24),  and  multiplying  the  four  densi¬ 
ties  together,  the  likelihood  function  (3.25)  is  formed. 


-1  4_ 

2  Y  (T.  -  U  +  DC./V) 

2  i  .  Cl  1  1 

1=1 


(3.25) 


In  (3.25)  C 4  has  teen  used  for  notational  :  onvenience, 
and  is  defined  to  be  zero.  Then  tae  log-like liaood  function 
is  formed  by  taking  natural  logs  of  (3.25),  yielding  (3.26). 


-2  In  :2  7T  )  -  4  In  !3)  - 


I  !V- °VV:'  (3-2«) 


Since  the  values  of  the  Ci  are  to  be  direction  cosines, 
the  usual  direction  cosines  constraint  must  be  added  to  LI 
to  form  the  Lagrangiar  function  L2  given  in  (3.27). 


1 


(3.27) 


=  l,  -  y  j  r>  fz:>  - 

-  A  Ln  1  J 

The  objective  is  to  maximize  L2  over  tie  possible 
choices  of  Cl,  C 2,  C3,  0,  S  and  lambda.  Ignoring  S  ior  the 
moment,  maximization  of  L2  can  be  acheived  by  minimization 
of  1  given  in  (3.28). 

14  •>  r  3  -  i 

“  '2  5-  I  1Ti  -  0  *  DCi'vr  'X  I,  IC?  "  1  (3.23) 

1=1  [l=l 

Take  partial  derivatives  of  L  to  get  (3.29)  i  nd  (3.30). 


T  (Ti  ‘  u  +  °VV)  '  -  \Ci  (3 .  25) 


ra  r  r  4  r,  3  1 

-  2  Vi.  -  40  .  -2-  V  c. 

Su  [in 1  v  hi  "J 


(3. 30) 


Eguate  (3.30)  tc  zero  and  solve  for  0  to  obtain  (3.3  1), 
the  maximum  likelihood  estimate  for  the  time  at  the  c  hy.dro- 
Fhone . 


/  4  n  3 

U  -  f  Jt.  *  -f-  T  c. 
V  ] = i  J  v  :*=i  ' 


(3-31) 


Eguate  the  three  equations  of  (3.29)  to  zero,  and 
multiply  each  equation  by  Ci  respectively,  to  ob;  ain  (3.32). 


X=i  ■  4-  (v. 


L'C  .  + 


4t) 


3  (3.32) 


Add  the  three  equations  of  (3.32)  together,  and  use  the 
direction  cosine  constraint  to  obtain  (3.33). 


X  - 


y  t.c.  -  u  y c.  +  — 

f->  l  1  f-  l  V 

1=1  1=1 


(3.  33) 


Equation  (3.34)  results  when  (3.29)  is  eguata  d  tc  zero. 


c. 

1 


T .  -  U 

1 


i  =  1,2,3 


(3.  34) 


Substitute  (-5.33)  into  (3.34)  to  obtain  (3.35),  the 
maximum  likelihood  estimate  of  the  direction  cosines,  where 
0  is  the  estimate  calculated  by  (3.31). 


T.  -  U 

ci  =  - J - - -  (3.35) 

y  T  .C  -  u  V  c  . 

Jr*,  j  ]  i 

j-1  :  =  1 

As  the  reader  will  perhaps  have  noticed,  tie  equations 
of  (3.35)  define  each  of  the  unknowns  Ci  in  terms  of  all 
three  unknowns.  Such  a  structure  suggests  that  (3.35)  can 
he  used  as  an  iteration  function.  That  is,  i:  reasonable 
initial  values  are  used  for  the  three  unknowns  in  the  right 
hand  side  of  (3.35)  then  new  values  are  produos  d.  Repeat 
the  process  using  the  new  values  until  the  answer  stabilises 
within  acceptable  tolerances.  Although  convergence  to  the 

correct  solution  is  not  guaranteed,  the  method  has  never 
failed  for  the  equations  of  (3.35).  Unlike  the  L.S.  method, 
(3.35)  dees  not  have  any  known  closed  form  solution. 

Returning  to  the  standard  deviation  S,  take  the  partial 
derivative  of  (3.27)  with  respect  to  S  to  get  (3.  36)  . 


3l2 

ds 


(3. 36) 


Multiply  (3.36)  by  S3  and  solve  for  S2 
maximum  likelihood  estimate  of  the  variance. 


to  get  the 
given  by 


(3.37). 

2 


r 


y  rr .  -  u> 

i=i  1 


(3. 37) 
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I.  THE  MAXIMUM  LIKE1IHOOD  PLAHAE  METHOD 

In  summary,  the  third  method  for  estimation  of  an 
apparent  position  is  as  follows: 

1.  let  U  =  T 4  initially; 

2.  use  the  L.S.  method  to  calculate  the  initial  values 
for  Ci,  (i=1,2,3),  using  (Se321); 

3.  use  U  and  Ci  (i=1,2,3)  in  the  right  hand  side  of 
{3.35)  to  obtain  new  estimates  for  Ci  (i  =  1,2,3); 

4.  recalculate  U ,  using  (3.31); 

5.  reiterate  steps  3  and  4  until  tne  values  Ci  (i=1,2,3) 
and  U  converge  within  acceptable  tolerances; 

6.  calculate  S2,  using  (3.37); 

7.  calculate  the  estimated  apparent  position  relative  to 
the  c-phone,  using  (3.33)  ; 


8. 


=  V  u  c. 


=  V  u 


■3  (3.33) 


lastly  translate  this  solution  and  its  c>  rresponding 
time  estimate  to  a  solution  and  time  relative  to  the 
acoustic  center,  using  (3.39),  where  E  and  Rc  are  as 
defined  by  (1.4)  in  Chapter  I. 


X  =  X  +  D/2  V  =  Y  +  D/2 

c  c 

T  =  U  R  /  R 


=  Z  +  D/2 
c 


(3.39) 


This  method  shall  be  hereafter  referred  to  as  the 
’Maximum  Likelihood  Planar  Method’,  or  M.L.P.  for  short. 
Originally  the  hopes  for  this  method  were  rather  high,  espe¬ 
cially  since  it  was  the  first  method  to  produce  a  variance 
estimate.  However  subsequent  experience  with  the  method 
indicates  that  it  probably  suffers  significantly  from  at 
least  two  factors: 

1.  the  planar  wavefront  assumption  probably  builds  in  a 
position  bias  as  in  the  case  of  the  L.S.  aethod;  and 


45 


2.  the  variance  estimate  is  inflated  since  part  of  tne 
noise  being  measured  is  due  to  the  inadequacy  of  the 
planar  assumption. 

J.  COMPUTATIONAL  EXAMPLE 

For  a  quick  comparison,  tne  three  aetnois  leveloped  in 
this  chapter  are  now  applied  to  the  example  which  was  used 
in  Chapter  II.  The  EXACT  results  calculated  previously  are 
included  in  Table  III  for  comparison. 


Method 

L.  S. 

L  .  S.  C . 

M.  L.  P. 

EXACT 


TABLE  III 

Single  Example  Comparison  of 
Planar  Wavefront  Methods 


Transit  Time 
7  7 secs) 

0.67529319 

0. 67527149 

0. 67529007 

0. 67527043 


Elev.  Angle 

niesir — 

15.22793 
15. 26372 
15. 10437 
15.27  002 


L  .  S.  C . 

a.  l. p. 

EXACT 


Apparent  Position  Estimate 
(  1003.  309  ,  3019.751  ,  866.  1!  6  ) 

(  1006.089  ,  3018.266  ,  868.235  ) 

(  997.722  ,  3023.685  ,  859.355  ) 

{  1005.966  ,  3017.899  ,  863.555  ) 


As  can  be  seen  easily,  the  M.L.  ?.  method  fares  rathe 
poorly  in  all  regards,  even  worse  t.ian  the  L.S.  method 
This  will  he  confirmed  by  the  evaluations  made  in  Chapter  V 
Also  worthy  of  note  is  tne  apparent  tendency  of  the  L.S.C 


IV.  A  SPHERICAL  HAVER FRONT  MODEL 

A.  TEE  SPHERICAL  WAVEFRONT  ASSUMPTION 

All  the  models  developed  in  Chapter  III  ire  limited 
primarily  by  the  assumption  that  the  wavefront  is  planar 
upon  arrival  at  the  hydrophone  array.  In  this  chapter  a 
model  is  developed  under  the  assumption  that  the  wavefront 
is  spherical  upon  arrival  at  the  array.  If  the  sound 
velocity  prefile  were  constant  with  depth  then  the  spherical 
model  would  be  exact.  This  is  of  course  not  the  case,  but 
it  is  suspected  that  the  wavefront  is  better  modelled  as  a 
sphere  than  as  a  plane  because  that  small  piece  of  the  wave- 
front  which  passes  through  the  30  foot  cute  spanned  by  the 
array  is  locally  spherical.  That  is  because  e/ ery  part  of 
that  piece  travelled  through  approximately  the  same  regions 
of  water,  experiencing  the  same  general  raytendng  patterns. 

The  spherical  assumption  is  accurate  if  and  only  if  the 
speed  of  sound  is  constant  over  the  ray  path,  and  conseg- 
uentiy  the  original  four  spherical  equations  apply  once 
again.  They  were: 

(  X  -  D/2  )2  +  (  Y  +  D/2  )2  -t-  (  Z  +  D/2  )2  =  V2T2 

x  (4.1j 

(  X  +  D/2  )“  +  (  Y  -  D/2  )2  +  (  Z  +  D/2  )2  =  y2T2 

y 

(  X  +  D/2  )2  +  (  Y  +  d/2  )2  +  (  Z  -  D/2  )2  =  Y2?2 

(  X  +  D/2  ) 2  +  (  Y  +  D/2  )  2  +  (  Z  +  D/2  ) 2  =  V2?2 

c 

It  has  been  stressed  previously  that  there  is  no  exact 
solution  (X,Y,Z)  satisfying  ail  four  eguations  (1.1).  That 
is  because  the  time  values  on  the  right  hand  side  correspond 
to  the  reality  of  a  variable  velocity  profile.  However,  if 
tne  spherical  wavefront  assumption  is  to  be  accurate,  then  a 
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constant  velocity  is  the  assumed  case,  and  any  ir.nacur  acies 
an  the  tame  values  are  regarded  as  due  to  timing  errors 
only.  Therefore,  under  the  spnerical  assumption,  if  the 
true  time  values  Ji  (i=1,2,3,4)  were  known  and  substituted 
into  (4.1),  an  exact  solution  to  the  overdeter  mined  system 
would  be  realized.  In  that  case  a  solution  to  ar.y  three  of 
the  equations  would  also  te  equal  to  taat  una  gue  exact  solu¬ 
tion.  In  particular,  the  UAVY  solution  of  Chapter  II  would 
be  the  true  solution.  Tnarefore  an  terms  of  tne  coordinate 
system  referenced  to  the  c  hydrophone,  X  would  te  given  by 
(4.  2)  . 


However,  X  is  also  given  by  (4.3)  ,  whers  Cl  is  the 
direction  cosine  along  the  X  axis  of  the  vector  from  the  c 
hydrophone  to  the  sound  source. 


(4.3) 


The  time  value  of  interest  at  the  moment  is  14,  the  time 
at  the  c  hydrophone.  Thererore  let  J  =  U4  for  clarity  of 
presentation,  and  then  equate  the  expressions  in  (4.2)  and 
(4.3)  in  order  to  solve  for  U.  If  this  same  lugic  is  also 
applied  to  the  similar  expressions  for  the  distances  to  the 
y  and  z  hydrophones,  the  results  are  (4.4). 


a  1.2,3  (4.4) 


The  expression  (4.4)  will  be  useful  in  the  development 
of  the  model  of  this  chapter. 


B.  LEAST  SCARES  MODELS 

A  direct  approach  might  be  to  apply  the  ljast  squared 
error  technique  to  the  spherical  equations,  an  a  manner 


paralleling  that  used  on  the  four  planar  equations  of  the 
L.5.  model  in  Chapter  III.  However  the  formulae  and  equa¬ 
tions  that  result  are  exceedingly  complex/  involve  fourth 
degree  powers  of  the  data  values  li,  and  have  thus  far 
defied  ail  solution  attempts.  Therefore  this  ida a  was  aban¬ 
doned  in  favor  of  the  maximum  likelihood  approach  winch 
follows. 


C.  MAXIMUM  LIKELIHOCE  COMPUTATIONS 

As  in  the  M.L.P.  model,  Gaussian  noise  is  assumed  for 
the  time  data  values.  Therefore 

Ti  =  Ui  +  Ei  i  =  1,2,3 

where  the  Ei  are  independent  identically  distributed  normal 
random  variables  with  mean  zero  and  variance  S2.  Tne 
density  of  each  Ti  is  therefore  (4.5)  . 

i  =  1,2, 3, 4  (a .5) 

To  form  the  likelihood  function,  multiply  the  four 
densities  together,  to  obtain  (4.6). 


Form  the  log  likelihood  function  ny  taking  nitural  loga¬ 
rithms  of  (4.6)  to  obtain  LI  of  (4.7)  . 


.  (T  .  )  = 

1  l 


n/-TT  s 


exo 


-1  ,  2 
- T  <  T.  -  U. 

-2  i  i 


-2  Ini 


-7T) 


(3) 


(4.7) 


In  order  to  maximize  Li  with  respect  to  Cl,  C2,  C3  and 
U,  it  is  sufficient  to  minimize  L2  of  (4.3),  where  a  substi¬ 
tution  for  Ui  has  been  made  using  (4.4). 
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(4.3) 


(T,-TJ) 


+  y 

1^1 


Ti  -  V'“‘  -  2. 


JC./V)  -  (l/v) 


Now  add  the  usual  direction  cosines  constraint  ar.i  born 
the  Lagrangian  function  1  of  (4.9). 

r  3 


-2  'A 


1=1 1 


-  1 


(4.9; 


For  notational  convenience,  let  Ki  be  given  by  (4.10). 


2  T 

-  U  -  (2  D  U  C.  /  V)  +  (D/V)“ 


.,2,3  (4.10) 


Take  the  partial  derivative  of  L  with  respect  to  Ci  to 
get  (4.11). 


"d  L 

3cx 


'2  ( T .  -  K.) 

1  i  -  2  D  U 


-2  K. 

l 


- 2  X  ci 


i  =  1,2,3 


(4.  11) 


Simplify  (4.11)  and  equate  to  zero  to  yield  (4.12). 

i  =  1,2,3  (4*  12> 


D  U  /  Ti 


C.  =  _i _ I 

1  X--  vy*~ 


Multiply  the  three  equation  in  (4.12)  by  Ci  (i=1,2,3) 
respectively,  and  add  them  togetner.  Thar,  use  the 
constraint  on  the  sum  of  squares  of  the  direction  cosines  to 
onta  in  (4.13). 


(4.  13) 


Substitute  (4.13)  into  (4.12)  to  get  the  max.  mum  likeli¬ 
hood  estimate  for  the  direction  cosines,  as  in  (4.14). 


5  1 


L  with  r;  spect  to  'J 


Now  take  the  partial  derivative  of 


(U 


T.)  (4-15) 


Equate  (4.15)  tc  zero  and  solve  for  LI  to  obtain  the 
maximum  likelihood  estimate  of  the  time  value  (4.16)  at  the 
c  hydrophone. 


u  = 


(4.  16) 


Finally  take  the  partial  aerivative  of  LI  in  (4.7)  with 
respect  to  £.  Equate  it  to  zero  and  solve  to  get  (4.17), 
the  maximum  likelihood  estimate  of  tne  variance. 


-> 


(T. 


(4.  17) 


D.  SOLUTION  BY  A  MODIFICATION  OF  NEWTON’S  METHOD 

The  solutions  given  by  equations  (4.14),  (4.16),  and 
(4.  17)  once  again  form  a  set  of  equations  which  would  seen 
to  be  solvable  by  natural  iteration  as  in  the  M...P.  method. 
Unfortunately  this  time  the  technique  fails  to  converge.  A 
mathematical  tool  is  needed  which  is  stronger  :  nan  natural 
iteration.  What  is  used  is  a  modiried  four  dimensional 
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version  of  Newton's  Method  [fief.  5  p.47]  to  sei  rch  for  the 
roots  of  d  set  of  four  equations. 

The  objective  is  to  determine  the  values  Cl,  C2,  C3  and 
U  which  satisfy  (4.14)  and  (4.16).  for  further  notational 
convenience  define  the  values  M  and  N  as  in  (4.  IS)  and 
(4.19)  . 


(4.  13) 


(4.  19) 


Given  those  definitions  of  M  and  N ,  then  tie  equations 
whose  roots  are  desired  can  be  simplified  to  (4.20). 


Now  define  error  functions  as  in  (4.21).  Ihise  evaluate 
the  amount  of  error  ir.  each  of  the  equations  (4.20)  for  any 
set  of  values  for  (C  1  ,C2  ,  C2  ,  0)  . 

(4.21) 

“*  * 

Let  G  he  the  four  di  me  ns  lor.ai  coil mn  vector 
(d  1  /  9  2,  <]3,  g  4)  '  .  Finally  1  *_•  t  d?  re  t..e  matrix  of  partial 
derivatives  of  G,  as  in  (4.  22)  . 

Newton's  Method  in  four  dimensions  says  t ha t  if  X  (n)  is 
a  four  dimensional  vector  noldiny  t:.o  current  approximate 
roots  Cl,C2 ,C3  and  U,  then  X(n+1)  will  ce  an  improved 
answer,  where  X  in*  1)  is  ji  v  en  by  (4.23). 
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Of  course  it  is  necessary  to  calculate  the  derivatives 
held  in  the  matrix  GP  in  order  to  use  (4.23).  Ti  ose  deriva¬ 
tives  are  given  in  Appendix  B, 

Unfortunately  when  multidimensional  versions  of  Newton's 
Method  are  applied,  there  is  often  a  tendency  for  the  method 
to  converge  slowly,  or  even  diverge.  This  is  tecause  it 
tends  to  overshoot  the  best  answer  for  each  iteration.  To 
alleviate  this  problem,  a  modification  is  made  to  the 
method.  At  the  end  of  each  Newton  iteration,  prior  to 
proceeding  with  the  next  iteration,  a  Golden  Section  Search 
[Ref.  6  ]  is  performed  to  find  the  best  possible  answer  in 
tne  direction  of  the  new  iterative  solution.  Specifically 
the  line  in  four  dimensional  space  from  X(n-1)  of  the 
previous  iteration  to  X(n)  of  the  present  i;  eration  is 
searched  for  the  best  answer.  The  definition  ot  the  'best' 
answer  is  that  point  along  the  search  line  which  minimizes 
the  sum  cf  squared  error  functions,  namely  (4.24|  . 


1  A 


( ».2«) 


The  current  iterative  solution  X(n)  is  thei  given  the 
value  of  the  minimizing  point  resulting  from  tne  Golden 
Section  Search.  Then  the  next  iteration  of  hew:  on' s  Method 
is  performed/  along  with  another  Golden  Section  Search  to 
find  the  next  iterative  solution  X(n+1). 

E.  THE  MAXIMUM  LIKELIHOOD  METHOD 

In  summary,  the  fourth  and  final  alternative  method  to 
estimate  apparent  positions  is  as  follows: 

1 .  let  U  =  T4  initially ; 

2.  use  the  L.S.  method  (Se321)  to  initially  estimate  the 
values  of  C  j  ( j=  1 , 2 ,  3 )  ; 

3.  set  X  ( 1 )  =  (C1,C2,  C3,  0)  ; 

4.  initialize  the  Newton  Iteration  counter  :  I  =  1  ; 

5.  calculate  the  values  Ki,  a  and  N  in  aco  rdance  with 
eguations  (4.10),  (4.13)  and  (4.19); 

6.  calculate  the  error  function  vector  G(X(I)),  using 
(4.20)  and  (4.21)  ; 

7.  calculate  the  derivative  matrix  GP(X(I)),  using  the 
results  in  Appendix  B; 

8.  invert  GP  (X(I))  ; 

9.  calculate  the  new  Newton  estimate  X(I+1)  from  (4.23)  ; 

10.  perform  a  Gclden  Section  Search  along  the  line 
hetween  X(I)  and  X(I+1)  to  find  the  point  which  mini¬ 
mizes  (4.  24)  ; 

11.  let  X  (I  +  1 )  be  egual  to  the  minimizing  point  found  in 
the  previous  step; 

12.  increase  the  Newton  iteration  counter  :  I  =  1+1  ; 

13.  reiterate  steps  5  through  12  until  the  values 

X{I)  =  (C1,C2,C3,U) 

converge  within  acceptable  tolerances; 

14.  calcualte  the  estimate  of  S2  using  (4.17). 
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This  Method  shall  hereafter  be  rtf«rrt-i  to  as  tne 
•  .lax  imam  Lixelihood  Spherica*  Method',  or  .1.1.5.  :or  short. 
As  will  re  seen  from  the  evaluations  made  in  o.d  ter  7,  tn*-. 
M.L.S.  method  is  apparently  the  only  alternative  to  consis¬ 
tently  rival  the  performance  of  the  currently  used  !I V  i _ a 

method.  It  has  the  additional  advantage  tnat  it  estimates 
the  error  present  in  the  time  data  values. 

F.  A  CGHPOTATIONAL  EXAMPLE 

This  latest  method,  I.L.S.,  is  now  applied  to  the  same 
example  considered  in  Chapters  II  and  III.  ?  or  a  yuick 
comparison.  Table  17  lists  the  results  ox  using  all  the 
methods.  The  error  vector  lengths  are  tne  distaaces  of  each 
position  estimate  from  the  EXACT  apparent  position. 

Since  this  is  only  one  example,  this  table  is  not 
presented  for  the  purpose  of  any  broad  conclusioas.  However 
it  is  of  interest  to  note  that  in  this  example 

1.  the  M.L.S.  method  outperforms  all  others,  including 
the  NAV  Y_  A  method;  and 

2.  in  many  ways  the  original  NAVY  metnod  outperforms  the 
adjusted  NA V Y_A  method. 

G.  VARIANCE  ESTIMATION 

In  the  computational  example,  the  time  data  values  used 
were  exact  since  the  methods  of  Appendix  A  coi  Id  be  used 
with  the  known  linear  velocity  profile.  Therefore  the 
appropriate  value  for  variance  in  the  time  valies  would  b<= 
zero.  For  this  error  free  example,  the  estimates  shown  in 
Table  V  were  obtained  for  the  standard  deviations  of  timing 
noise,  using  the  two  maximum  likelihood  methods. 

In  the  case  of  this  one  example,  the  spherical  assump¬ 
tion  is  apparently  an  improvement  over  the  planar,  since  the 
M.L.S.  estimate  of  error  is  only  4%  of  the  M.L.P.  estimate. 
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I  ABLE  17 

Single  Example  Comparison  of  All  Methods 


Method 

Transit  Time 

1  T§ecs) 

Elev.  Angle 

FTIegfr 

L  engt  h 
Error  ; 

of 

Vector 

NAVY 

0.  67526969 

1  5.  26459 

0.412  8 

NAVY_A 

0. 67527027 

1  5.  26461 

0.4840 

L.  S. 

0.  67529319 

15.22798 

3.730 

L  .  S.  C . 

0. 67527149 

15. 26372 

0.522 

M.  L.  ?. 

0. 67529007 

15.  10437 

13.64) 

M  .  L .  5  . 

0. 67527C49 

1 5.26677 

0.4117 

EXACT 

0.  675 2  7043 

15.  27002 

A  p parent 

Position 

2 

stimate 

NAVY 

(  1005.957 

,  3017.874 

/ 

868. 11 3 

) 

N AVY_A 

(  1  006.  087 

,  3018.259 

t 

868.25  5 

) 

L.  S. 

{  1  003.  809 

,  3019.751 

t 

86  6.  12  6 

) 

L.S.  C. 

{  1  006.  089 

,  3018.266 

t 

868.235 

) 

a. i. p. 

(  997.722 

,  3023.635 

/ 

859.35  5 

) 

M. 1. S. 

{  1006.195 

,  3018.190 

t 

868.375 

) 

EXACT 

(  1005.  966 

,  3017.899 

t 

868.  55  5 

)' 

That  is  reyarded  as  an  improvement  because  the  correct 
answer  is  zero.  The  higher  M. L. P.  estimate  is  indicative  of 
the  inflation  due  to  the  planar  assumption. 

The  tracking  range  which  provided  data  for  this  study 
records  time  values  to  seven  decimal  places.  Tierefore  the 
standard  deviation  estimated  by  the  M.L.3.  method  is  of 
particular  interest,  since  it  indicates  errors  in  the 
seventh  decimal  place  even  when  there  is  no  error  present. 
Since  the  data  was  actually  error  free  in  this  ecampie,  the 


TABLE  V 

1  Maximum  Likelihood  Error  Estimates 

for  the  Example 

Problem 

Model  Method 

1st.  Std.  Devil  tion 

Planar  M.L.P. 

5.517  E-6  secs. 

Spherical  M.L.S. 

2.191  E-7  se:  s . 

estimate  is  a  measure  of  the  variation  induced  b{  the  spher¬ 
ical  wavefront  assumption.  A  broader  discussion  of  error 
estimation  will  be  undertaken  in  Chapter  VI. 
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V.  IliLOATICNS  OP  MODELS 


A.  GENERAL 

There  are  many  sources  of  errors  in  the  ovsrall  hydro- 
phonic  tracking  problem.  These  include,  among  others: 

1.  errors  in  estimation  of  the  sound  velocity  profile; 

2.  inho mogeneity  of  the  velocity  profile  over  time  and 
horizontal  displacement;  and 

3.  possible  errors  in  measuring  the  positioi s,  and  the 
angles  of  tilt  and  rotation  for  the  hydrophonic 
arrays. 

This  study  focuses  on  those  errors  which  occur  during 
the  computations  preceding  the  ray  tracing  procedure.  To 
evaluate  the  performance  of  the  methods  developed  in 
Chapters  II,  III  and  IV,  it  is  necessary  to  control  strictly 
the  ray  tracing  procedure.  Only  in  that  way  can  the  differ¬ 
ences  found  between  methods  be  attributed  to  the  differences 
between  models,  and  not  to  any  source  of  error  oj tside  those 
met  hods . 

It  jas  originally  hoped  that  the  various  methods  mignt 
be  compared  by  applying  them  to  real  tracking  data.  However 
it  was  found  that  the  overall  tracking  problem  iad  too  many 
iurjfc  sources  of  error  to  allow  the  methods  to  demonstrate 
a;..-  Ii  : :  <=r r.ces.  Therefore  the  methods  were  con  pared  under 
a  mol*,  t.  ij.ntly  controlled  simulated  environment. 

E.  SITUATION  SCENARIO 

Two  ;:;n:er  t  simulations  are  used  to  compare  the  six 
dr::>  :*  t  i*  tnods.  Eoth  use  a  basic  scenario  similar  to 
that  c  t  *■...*  compu  ta  tional  example  explored  in  Chapters  II, 
III  an:  IV.  That  example  assumed  tnat: 
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1.  the  velocity  versus  depth  relationship  is  linear,  and 
given  exactly  ty: 

V  =  4840.7  +  0.03314  •  DEPTi  ; 

2.  the  acoustic  center  of  each  array  is  at  a  depth  of 
1300  feet; 

3.  the  hydrophone  arrays  are  ail  level,  and  their  X,  Y 
and  Z  arms  are  parallel  to  the  respective  coordinate 
axes  of  the  tracking  range. 

Under  these  circumstances  the  methods  se:  forth  in 
Appendix  A  can  be  used  to  compute  the  exact  values  for  the 
hydrophone  times,  ray  transit  time  and  elevation  an^le  for  a 
sound  wave  emanating  from  a  source  at  any  specified  loca¬ 
tion.  Therefore  when  the  methods  are  applied  to  those  exact 
times,  the  resulting  estimated  positions  can  be  compared  to 
the  known  true  position. 

C.  SIMULATED  ERROR  FOR  TIMING  DATA 

The  models  developed  in  this  study  were  designed  to 
improve  position  estimation,  especially  in  the  presence  of 
errors  in  the  timing  data.  Lacking  any  better  m3  del  at  this 
time  for  timing  errors,  the  simulated  environment  includes 
an  assumption  of  Gaussian  errors  for  the  hydrophone  times. 
Therefore  realistic  timing  data  can  be  simulate!  by  adding 
to  each  exact  time  value  a  random  guantity  of  normally 
distributed  error.  The  mean  of  the  error  is  assumed  to  be 
zero.  The  variance  was  estimated  from  real  tricking  data, 
using  the  variance  estimating  property  of  the  M.L.S.  model. 
The  data  from  one  tracking  run  was  used,  involving  six 
hydrophone  arrays  and  733  position  estimates  (sea  Chapter  VI 
for  data  selection  details).  Each  position  from  the 
tracking  run  produces  one  estimate  of  the  variance.  The 
variance  Vdxue  chosen  for  use  in  the  simulations  was  the 
median  of  the  733  variance  estimates  produced  b f  the  M.L.S. 
method.  That  value  was 


S2  =  9.  1204  2-12  secs2  . 

That  is  the  same  as  a  standarad  deviation  of 

S  =  3. 02  2-6  secs  . 

Each  of  the  two  types  of  simulations  was  r  ui  four  sepa¬ 
rate  times.  Each  run  was  done  with  a  different  specified 
variation.  Those  four  distinct  error  conditions  are  delin¬ 
eated  in  Table  VI  . 


D.  SINGLE  ARRAY  SIMULATION 

In  the  first  simulation,  the  intent  was  to  compare  the 
methods  pairwise,  sc  as  to  determine  which  me;  hod  is  more 
likely  to  produce  the  more  accurate  estimate  of  a  given 
sound  source  position.  One  thousand  positions  were  chosen 
at  random.  Each  position  was  3000  feet  from  the  array.  The 
positions  were  uniformly  spread  over  the  surface  of  a  sphere 
of  radius  3000  feet,  centered  at  the  array,  truncated  above 
by  the  water  surface  (depth  0}  and  below  by  the  depth  of  the 
array  (1300  feet).  The  methods  set  forth  in  Appendix  C  were 
used  to  assure  that  the  random  positions  selected  were 
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uniformly  distributed  over  the  surface  area  of  tne  truncated 
hemisphere. 

Each  of  the  1000  randomly  chosen  source  pusitions  was 
then  processed  as  follows; 

1.  calculate  the  exact  hydrophone  times,  using  the 
methods  of  Appendix  A; 

2.  add  to  each  of  the  four  exact  times  a  random  value  of 
error  at  the  specified  level  (zero,  low,  medium  or 
high) ; 

3.  apply  each  of  the  six  metnods  to  the  hydrophone 
times,  generating  six  different  apparent  positions; 

4.  apply  the  ray  tracing  procedure  to  each  of  the  six 
positions,  using  layers  that  are  25  feet  thick,  and 
utilizing  the  known  linear  velocity  profile,  thereby 
producing  six  different  estimates  of  the  sound  source 
lcca  tion ; 

5.  compare  the  six  different  estimates  pairwise  to  see 
which  method  in  each  pair  produced  the  estimate 
closest  to  the  true  sound  source  location. 

The  comparison  being  made  is  that  one  method  is  consid¬ 
ered  preferrable  to  the  other  if  it  more  frequently  produces 
the  more  accurate  estimate. 

The  layer  thickness  of  25  feet  was  selected  order  to 
simulate  the  actual  procedure  at  the  tracking  range  which 
provided  data  for  this  study.  however,  as  discussed  in 
Chapter  I,  when  the  velocity  profile  is  smooth  and  known 
exactly,  the  process  is  very  robust  with  respect  to  the 
thickness  used.  The  thickness  values  1,  10,  and  25  were  all 
tried,  with  virtually  no  changes  in  any  of  the  comparisons 
between  methods. 
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E.  SINGLE  AfiRAY  SIMULATION  RESULTS 

Tables  VII  and  VIII  contain  the  results  of  the  single 
array  simulation.  Each  tabled  entry  represents  :  ne  fraction 
of  time  that  the  method  of  that  row  produced  a  better  esti¬ 
mate  than  the  method  of  that  column.  For  example,  in  Table 
VII,  the  L.S.  method  outperformed  the  M.L.S.  me:  hod  in  only 
5.35;  of  the  1300  trials  with  low  error  values. 

In  a  one  tailed  test  that  one  method  is  setter  than 
another,  tnese  binomial  proportions  are  significant  at  the 
0.05  level  if  they  exceed  0.526.  Symmetrically,  one  method 
is  significantly  worse  than  another  at  the  0.05  level  if  the 
proportion  is  less  than  3.474.  For  the  0.01  level  the 
corresponding  critical  values  are  0.537  and  0.463 
respectively. 

The  results  indicate  that: 

1.  as  previously  claimed,  the  NAVY_A  method  usually 
outperforms  the  unadjusted  NAVY  method;  of  particular 
interest  is  the  case  of  zero  error  vhi; h  actually 
compares  the  relative  ability  of  eacn  method  to 
produce  the  exact  answer  when  given  the  exact  times; 
in  those  cases  the  N  AV Y_A  method  does  extremely  well 
against  the  NAVY  method; 

2.  under  all  error  conditions  the  most  successful 
performer  is  the  M.L.S.  method,  since  it  always  has  a 
favorable  (greater  than  0.5)  comparison  fraction 
against  all  other  methods;  the  M.L.S.  fractions  vary 
little  over  the  four  error  levels; 

3.  under  all  error  conditions,  the  worst  performer  is 
always  the  M.L.P.  method; 

4.  spherical  methods  consistently  outperform  planar 
methods;  and 

5.  increased  error  levels  tend  to  lessen  the  distinction 
between  methods;  in  the  zero  error  case  comparisons 
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TABLE  VII 

Single  Array  Simulation  Results 
for  Lower  Error  Levels 

ERROR  LEVEL  :  ZERO 

NAVY  N  AV Y_ A  L.S.  L.S.C.  M.L.P. 

M  T  9 

u  •  1_  m  • 

NAVY 

.011 

.950 

.676 

.987 

.  420 

N  A VY_A 

.98  9 

.95  0 

.728 

.987 

.434 

L.  S. 

.050 

.  050 

.050 

.937 

.  057 

L.  S.C. 

.32  4 

.  272 

.950 

.987 

.  346 

M.  L. ?. 

.01  3 

.  013 

.01  3 

.013 

.013 

A. L. S. 

.560 

.  566 

.943 

.6  54 

.987 

NAVY 

ERROR 

N  A VY_A  L 

LEVEL  : 

.S.  L 

LON 

.S.C. 

M.L.P. 

a. l.s. 

NAVY 

.  165 

.952 

.648 

.987 

.396 

N  A  VY_A 

.835 

.952 

.693 

.987 

.  405 

L  .  S  • 

.048 

.  048 

.043 

.  987 

.058 

L.  S.C. 

.352 

.  307 

.952 

.987 

.  369 

M.L. P. 

.01  3 

.  013 

.013 

.013 

.013 

a.i.s. 

.604 

.  595 

.  942 

.631 

.987 

are  in  the  interval  (0.01,0.99),  while  in  the  high 
error  case  that  interval  is  narrowed  considerably  to 
(0.  37,0.63)  . 

F.  DOUBLE  ARR  A  Y  SIflOIATION 

In  the  second  simulation,  the  intent  again  was  to 
compare  the  methods  pairwise,  this  time  determining  which 
method  is  mere  likely  to  produce  positions  which  agree  more 
closely  in  the  two  array  crossover  problem.  This  is  not  the 
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TABLE  VIII 

1 

Single  Array 
for  High 

Simulation  Results 
er  Error  Levels 

ERE  OF. 

LEVEL  : 

MEDIUM 

NAVY 

N  AVY_ A 

L  .  S.  L 

.  s .  c . 

M.L.?. 

M .  L .  3  . 

NAVY 

.  464 

.819 

.559 

.987 

.  435 

N  A  V Y_A 

.54  6 

.821 

.566 

.937 

.  437 

L.  S. 

.131 

.  179 

.181 

.971 

.  177 

L. S. C. 

.44  1 

.  434 

.819 

.937 

.  450 

M.  I. F. 

.013 

.  013 

.  029 

.0  13 

.013 

3.1. S . 

.56  5 

.  563 

.  323 

.550 

.937 

ERROR 

LEVEL  : 

HIGH 

NAVY 

N  A VY_A 

L.  S.  L 

.  S.  C. 

M.i..  P . 

M  .  L  .  S  . 

NAVY 

.  486 

.  544 

.439 

.562 

.  431 

N  A VY_ A 

.51  4 

.546 

.516 

.559 

.  431 

m  S  m 

.456 

.  454 

.457 

.557 

.  400 

L.S.C. 

.56  1 

.  484 

.543 

.  5  g2 

.  427 

M.l. P. 

.438 

.  441 

.  443 

.438 

.  373 

3.1.2. 

.56  9 

.  569 

.600 

.573 

.  627 

same  Question  as  that  addressed  by  the  single  array  simula¬ 
tion.  The  two  estimates  produced  by  any  one  m2  thod  may  be 
very  close  to  each  other,  and  yet  be  iar  away  from  the  true 
position. 

For  the  double  array  scenario  two  arrays  are  used,  sepa¬ 
rated  by  7500  feet,  both  at  depths  of  1300  feet..  The  arms 
of  each  array  are  parallel  to  the  corresponding  coordinate 
system  axes.  Once  again  1000  positions  were  randomly  gener¬ 
ated  in  a  uniform  manner,  this  time  over  a  3  dimensional  box 
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Figure  5.1  Double  Array  Siaulation  Configuration. 

running  crossways  between  the  two  arrays.  The  box  is  1300 
feet  deep,  5000  feet  long  and  1000  feet  across  (see  figure 

5.1)  . 

Each  of  the  1000  randomly  chosen  sound  source  positions 
in  the  two  array  simulation  were  processed  as  foLlows: 

1.  calculate  the  exact  hydropnone  times  for  the  first 
array  using  the  methods  set  forth  in  Appeidix  A; 
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2.  add  to  the  exact  times  some  random  error  at  the  spec¬ 
ified  level; 

3.  Repeat  steps  1  and  2  for  the  second  array,  usiny  a 
new  set  of  random  error  values; 

4.  apply  each  of  the  six  methods  to  the  time  values  from 

both  arrays,  producing  six  differen:  pairs  of 

apparent  positions; 

5.  apply  the  ray  tracing  procedure  to  bo; h  apparent 
positions  in  each  of  the  six  pairs,  utilizing  the 
known  linear  velocity  profile,  producing  six  pairs  of 
estimated  sound  source  positions; 

6.  for  each  of  the  six  pairs  of  positions,  ciiculate  the 
distance  between  the  two  positions  in  the  pair; 

7.  make  pairwise  comparisons  of  the  six  different 

distances,  to  see  which  method  in  each  pair  exhibits 
the  closest  agreement  between  its  two  position  esti¬ 
mates. 

This  time  the  comparison  being  made  is  tha;  one  method 
is  considered  better  than  anotner  if  the  positions  it 
produces  agree  more  closely  more  often  than  those  of  the 
other  method. 

G.  DCOBLE  ARRAY  SIMDIATION  RESULTS 


Tables  IX  and  X  contain  the  results  of  the  1 ouble  array 
simulation.  Each  tabled  entry  represents  the  fraction  of 
time  that  the  method  of  that  row  produced  a  pair  of  esti¬ 
mates  which  were  in  closer  agreement  than  the  estimates 
produced  by  the  method  of  that  column.  Significance 
criteria  for  these  proportions  are  the  same  as  for  the 
single  array  simulation. 

These  results  are  similar  to  those  of  the  single  array 
simulation,  because  thev  indicate  that; 


TABLE  IX 


Double  Array  Simulation  Results 
for  Lower  Error  Levels 


EE  EOF. 

LEVEL 

:  ZERO 

NAVY 

N  A  V  Y_A  L 

.5. 

L. S. C. 

M.  L. P. 

M.L.S 

NAVY 

.  366 

.989 

.726 

.  968 

.212 

N  A VY_ A 

.634 

.389 

.734 

.968 

.212 

Xi  •  0  • 

.01  1 

.  Oil 

.011 

.555 

.014 

L. S.C. 

.27  4 

.  266 

.989 

.968 

.  193 

M.L. P. 

.03  2 

.  032 

.  445 

.032 

.033 

M.L.  S. 

.78  8 

.  788 

.986 

.807 

.967 

ERROR 

LEVEL 

:  LON 

NAVY 

N  A VY_A  L 

.  S. 

i. .  S  .  c  . 

M.L.?. 

il .  L .  S 

NAVY 

.  448 

.989 

.661 

.952 

.310 

N  A VY_ A 

.55  2 

.983 

.684 

.952 

.  312 

L.  S. 

.01  1 

.011 

.011 

.560 

.014 

L  .  S.  C . 

.339 

.216 

.939 

.952 

.  254 

M.  L.  ?. 

.04  8 

.048 

.  440 

.048 

.038 

M.L.S. 

.690 

.  688 

.986 

.746 

.962 

the  H  AV  Y  _  A  method  outperforms  the  original  unadjusted 
NAVY  method,  although  the  difference  is  not  signifi¬ 
cant  in  the  higher  error  level  cases; 

the  most  successful  performer  is  consistently  the 
M.L.S.  method; 

spherical  methods  almost  always  outperform  planar 
methods ; 

increased  error  levels  tend  to  lessen  the  distinction 
between  methods. 


TABLE  X 

Double  Array  Simulation  Results 
for  Higher  Error  Levels 


ERROR 

LEVEL 

:  MEDIUM 

NAVY 

N  A  V  Y_  A  L 

.  S. 

L.S.C. 

M.L.?. 

M.I.S 

NAVY 

.  498 

.  84  1 

.513 

.308 

.  443 

N  A  VY_ A 

.53  2 

.  842 

.542 

.809 

.  440 

L-  S. 

.159 

.  153 

.159 

.525 

.  160 

L - S. C. 

.437 

.  458 

.84  1 

.808 

.  433 

M.L.  P. 

.192 

.  191 

.475 

.192 

.  144 

M.  1.  5. 

.557 

.  560 

.840 

.567 

.856 

ERROR 

LEVEL 

:  HIGH 

NAVY 

N  A VY_A  L 

-S. 

L.S.C. 

M.L.?. 

M.I.S 

NAVY 

.  487 

.556 

.535 

.479 

.  442 

N  A VY_ A 

.513 

.558 

.499 

.481 

.  440 

L.S. 

.444 

.  442 

.445 

.451 

.422 

L. S. C. 

.465 

.  501 

.555 

.480 

.438 

M.L. ?. 

.52  1 

.  519 

.54  9 

.520 

.  428 

M.I.S. 

.55  8 

.  560 

.578 

.562 

.  572 

There  are  however  sor.e  indications  from  these  results 
which  are  different  from  those  of  tne  single  array  simula¬ 
tion,  such  as: 

1.  the  worst  performer  was  consistently  the  L.S.  method, 
rather  than  the  M.L.  P.  method; 

2.  in  the  nigh  error  case  the  M.L.P.  method  ls  equal  to, 
or  even  marginally  better  than,  every  other  method 
except  the  M.I.S.  method; 


3.  the  perf ormance  of  the  M.L.  S.  method  i>  noticeably 
better  under  error  free  conditions. 

The  last  two  inferences  are  pernaps  the  most  inter¬ 
esting.  The  maximum  likelihood  approach  was  intended  to 
estimate  and  account  for  Oaussian  errors  in  the  timing  data 
values.  Hence  it  is  really  not  very  surprisii  g  that  the 
M.L.S.  method  does  well  with  error  prone  data.  3ut  it  is 
interesting  to  note  that  the  M.L.P.  method  also  does  well 
under  high  error  levels,  even  though  it  probably  suffers 
from  a  bias  due  to  the  planar  assumption. 

Cf  even  greater  interest  is  that  the  M.L.S.  method  seems 
to  Le  at  its  best  when  compared  to  otner  methods  under  error 
free  conditions.  This  result  was  unexpected,  and  indicates 
that  the  M.L.S.  method  not  only  handles  timing  errors  well, 
as  was  intended,  but  apparently  also  does  an  evei  better  gob 
of  approximating  the  elusive  exact  solution  to  the  original 
four  spherical  eguations  of  (4.1)  and  (2.1)  when  the  exact 
time  values  are  available. 


H.  LIMITATIONS  ON  I NTERPBETATION  OF  RESDLTS 

The  results  of  both  simulations  seem  to  imply  that  the 
M.L.S.  method  outperforms  the  currently  used  N A VY_A  method 
on  any  randomly  chosen  sound  source  position,  with  or 
without  timing  errors.  These  are  encouraging  results. 
Nevertheless  the  reader  is  cautioned  that  thesi  tests  are 
just  simulations,  and  like  all  simulations,  must  make 
assumptions  which  cannot  fully  reflect  the  reali;  y  of  actual 
hydrophonic  tracking  conditions.  The  most  important  assump¬ 
tions  made  for  these  simulations  are: 

1.  the  sound  velocity  profile  is  known  exactly,  and  is 
linear; 

2.  the  errors  in  the  timing  values  from  any  nydrophone 
are  normally  distributed  with  mean  zero,  and  are 
independent  of  the  noise  in  any  other  hydrophone;  and 
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the  parameters  of  the  simulation  were  fixed;  for 
example  the  single  array  simulation  used  a  fixed 
range  of  3000  feet,  and  both  simulations  used  arrays 
at  13C0  feet  of  depth,  with  fixed  oriental  ions  to  the 
range  coordinate  system. 

The  first  assumption  is  probably  of  little  consequence. 
It  not  only  greatly  facilitates  computations,  but  also  helps 
to  isolate  the  initial  angle  and  time  estimation  problem 
from  the  unrelated  errors  involved  in  the  velocity  profile 
estimation  procedure. 

The  second  assumption  is  somewhat  more  troublesome. 
Errors  may  not  be  Gaussian  at  ail,  or  if  they  are,  the  mean 
may  not  be  zero.  Unfortunately  each  positioi  estimation 
involved  only  four  equations,  and  therefore  did  not  allow 
for  estimation  of  more  than  four  paranenters.  Iierefore  the 
error  mean,  being  a  fifth  parameter,  could  not  be  estimated. 
Also  the  errors  of  any  one  hydrophone  may  very  well  not  be 
independent  of  the  errors  of  the  other  three  pha nes  on  that 
array.  fortunately  these  concerns  are  offset  somewhat  by 
the  results  of  the  double  array  simulation,  wherein  it  was 
found  that  the  K.L.S.  method  was  at  its  best  wnen  there  was 
no  noise  at  all. 

Also  the  error  type  and  level  may  depend  on  other 
factors,  such  as  the  target's  range,  elevation  and  azimuth 
angles  from  the  array.  Th is  highlights  the  concerns  of  the 
third  assumption.  There  is  considerable  room  here  for 
future  work  concerning  the  dependency  of  results  on  such 
complicating  factors. 

Lastly  it  should  be  pointed  out  that  the  simulations 
make  comparisons  only  on  the  binomial  basis  of  better  versus 
worse  in  1000  trials.  The  magnitudes  of  the  actual  differ¬ 
ences  are  ignored.  It  is  possible,  though  perhaps  unlikely, 
that  while  one  method  marginally  outperforms  a  second  method 
in  most  trials,  in  all  the  remaining  trials  the  :  irst  method 
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VI.  COHCIOSIONS  AND  HECOMdEN DATIONS 


A.  ESTIMATION  OF  TIMING  DATA  VARIATION 

One  of  the  primary  purposes  of  this  study  i  as  to  esti¬ 
mate  the  amount  of  variability  in  the  timing  data  being 
recorded  during  actual  tracking  runs.  This  problem  was 
addressed  by  the  M.L.P.  and  M.I.5.  maximum  likelihood 
models.  The  M.L.P.  model,  as  previously  discussed,  suffered 
from  a  bias  due  to  the  planar  assumption,  was  the  poorest 
estimator  of  positions  among  all  the  models,  and  produced  an 
inflated  variance  estimate.  Therefore  the  spherical  model 
M.L.S.  is  used  to  estimate  the  data  variance. 

The  variability  that  is  measured  by  the  M.L.S.  model  is 
made  up  of  three  components.  First  there  is  :  he  variance 
induced  by  the  spherical  assumption.  Then  there  is  the 
variance  caused  by  the  seven  decimal  accuracy  used  when 
recording  the  data.  Finally  there  are  the  errors  inherent 
in  the  physical  process,  due  to  such  factors  as  hydrophone 
variability  or  malfunction,  local  distortions  of  the  sound 
wave,  and  inexactness  of  the  water  column  which  estimates 
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Therefore  the  data  variance  can  oe  estimated  ty  nrst 
estimating  the  variability  induced  by  the  spherical  assump¬ 
tion,  and  then  subtracting  it  from  the  M.L.S.  estimate  of 
the  variation. 

The  M.L.S.  estimate  was  obtained  by  applying  the  M.l.S. 
method  to  the  data  from  a  tracking  run  at  the  Nar.ocse 
torpedo  tracking  range  on  May  6,  1980.  That  run  involved 
position  estimation  ty  several  different  hydrophone  arrays. 
The  run  made  several  thousand  position  estimates,  733  of 
which  were  at  depths  of  100  feet  or  more  and  involved 
targets  not  more  than  4700  feet  from  the  sensing  array.  The 
depth  limitation  was  imposed  to  avoid  the  excessive  compli¬ 
cations  caused  by  the  radical  changes  in  the  velocity 
profile  above  that  depth  (see  figure  2. 1) .  The  maximum 
range  limitation  imitates  the  data  validation  procedure  at 
the  tracking  range,  where  positions  farther  than  4700  feet 
from  the  array  are  discarded. 

The  tracking  data  yielded  the  following  range  of  esti¬ 
mates  for  the  standard  deviation  of  the  data  noise,  using 


M.L.S. 

model. 

MAXIMUM  VAIUE 

2.89 

E-5 

secs 

0.95  QUANTILE 

1.  18 

E-5 

secs 

MEDIAN  VAIUE 

3.  02 

E-6 

secs 

0.05  QUANTILE 

4.  1 1 

E-7 

secs 

MINIMUM  VAIUE 

3.  13 

E-8 

secs 

For  an 

overall  estimate  of 

the  noise. 

the 

was  used,  so  that: 


(  3.02  S-6  )  2 


The  variability  induced  by  the  spherical  assumption  was 
estimated  ty  applying  the  M.l.S.  procedure  to  perfectly 
noiseless  data  in  the  idealized  environment  of  the  single 
array  simulation  of  Chapter  7.  This  was  done  for  targets  at 
ranges  of  1500  to  450C  feet,  at  500  feet  increments,  with 


1000  randomely  chosen  targets  at  each  range 
are  collected  in  Table  XI  . 


The  results 


TABLE  XI 

Inflation  of  Error  Estimates  Induced 
by  the  Spherical  Model 


Standard  Deviation  Estimates  from  Exact  Tims  Values 
with  Known  Linear  Velocity  Profile 


P.  ANGZ 

MINIMUM 

Q  (-05) 

MEDIAN 

Q  (-95) 

MA{ I MUM 

1500 

2.122-  10 

2.84E-8 

2. 08E-7 

3. 13E-7 

3. 2 9  E-7 

2000 

2.12E-9 

5 . 88E-8 

2. 26E-7 

3. 14E-7 

3.3  0E-7 

2500 

3.90E-8 

1  .09 E-7 

2. 34E-7 

3.14E-7 

3.3  IE-7 

3000 

8.35S-8 

1 .50E-7 

2. 37E-7 

3. 15E-7 

3.3  4E-7 

3500 

1.21E-7 

1 . 59 E-7 

2.37E-7 

3. 13E-7 

3. 24E-7 

4000 

1 .47E-7 

1  .65E-7 

2.39E-7 

3. 16E-7 

3.2  5E-7 

4500 

1 .46  E-7 

1  .69 E-7 

2. 42E-7 

3.14E-7 

3. 25E-7 

Table  XI  shows  that  the  median  and  maximum  inflation 

values  are  reasonably  independent  of  target  range.  The 

minimum  values  vary  somewhat,  tut  only  for  nnge  values 

below  3000  feet.  This  represents  a  very  stable  situation 

overall.  Therefore  the  inflation  due  to  the  spherical 

assumption  is  estimated  by  the  median  value  at  a  range  of 

3000  feet,  namely: 

rr2  u  =  (  2.37  E-7  )  2  . 

t/  sph 

Combining  these  two  estimates,  the  variar.;e  estimated 
for  the  timing  data  is 


u  time 


ml  s 


u  sph 


(3.322-6)2  - 

(  3.31  1  E-6  )  2 


(  2.37  2-7  ) 


As  car.  be  seen,  the  error  induce!  by  the  spherical  model 
is  less  than  10/?  of  the  H.  L.S.  error  estimate.  Therefore 
when  it  is  accounted  for  by  subtraction  from  the  H. L.S. 
variation  estimate,  the  final  variance  estimite  changes 
little. 

The  estimated  value  indicates  a  standard  error  in  the 
6th  decimal  place.  With  a  typical  speed  of  sound  value  of 
4880  feet  per  second,  this  represents  a  position  differen¬ 
tial  of  about 


4380  •  3 . 01 1 E-6 


3.015  feet 


This  estimate  is  quite  low,  indicating  that  the  time  values 
being  recorded  are  sufficiently  accurate. 

There  is  considerable  opportunity  for  addi:  ional  work 
determining  the  relationship,  if  any,  between  the  time  vari¬ 
ance  and  other  factors  such  as  angles  of  elevation  and 
azimuth  of  the  target  from  the  array. 

There  is  also  the  problem  that  the  time  /ariatior.  is 
likely  tc  be  array  dependent.  For  example,  consider  figure 
6.1  ,  wherein  the  standard  error  estimates  from  the  actual 
tracking  run  are  plotted  versus  the  range  of  the  target. 
The  plot  does  not  indicate  that  there  is  any  simple  rela¬ 
tionship  between  range  and  error  level.  Howevjr  the  plot 
does  show  a  bunching  pattern.  When  the  error  estimates  are 
plotted  separately  for  each  array,  tnen  the  bunching  pattern 
becomes  clearly  associated  with  the  individual  arrays. 
Consider  figures  6.2  and  6.3,  where  the  separate  plots  have 
been  made  for  four  different  arrays.  It  is  still  not  clear 
from  these  plots  whether  the  principle  effect  is  due  to  the 
individual  arrays,  or  the  ranges  of  the  targets.  However 
some  level  of  array  dependency  seems  likely,  indicating  a 
need  for  additional  investigation. 


Figure  6.1  Error  Estimation  Versus  Range  of  Target. 

B.  CHOICE  OF  METHOD 

Clearly  all  indications  are  that  the  plani r  wavefront 
models,  L.S.  and  M.L.F.  are  not  candidates  for  use  as  posi¬ 
tion  estimators.  Furthermore  the  hybrid  model  L.S.C.  is  an 
interesting  improvement,  but  never  really  performs  well 
enough  compared  to  the  M.L. S.  and  NAVY  models. 

The  original  NAVY  model  is  usually  outperformed  by  the 
adjusted  NAVY_A  method.  However  the  differences  are  not 
always  significant. 

The  spherical  model  M.L. 5.,  on  the  other  haid,  consis¬ 
tently  outperformed  all  other  methods  during  the  simulated 
evaluations.  It  wculd  seem  that  M.L .  S.  is  tie  model  of 
choice.  It  does  the  best  job  of  handling  normally  distrib¬ 
uted  errors  in  the  data.  But  that  is  not  tie  strongest 


DATA  FROM  ARRAY  NO.  7 


DATA  FROM  ARRAY  NO.  14 


Figure  6.2  Error  Estimation  for  Arrays  7  and 
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STD.  DEV.  (10-5  SECS  )  STD. 


Figure  6.3  Error  Estimation  for  Arrays  56  and  57 
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A  sore  important. 


ana  surprising. 
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argument  for  its  use. 
argument  in  its  favor  is  that  when  the  exact,  error  free 
times  are  used,  the  h.L.S.  apparent  position  estimate  will 
usually  produce  the  most  accurate  estimate  of  the  true 
apparent  position.  This  is  the  desired  overall  result,  so 
that  the  sensitive  ray  tracing  procedure  will  be  affected 
minimally  ty  apparent  posit  ion  estimation  errors.. 

There  are  nevertheless  several  notes  of  caution  which 
should  be  considered  before  embracing  the  H.L.S.  method 
wholeheartedly.  The  first  caution  has  been  stressed  before, 
namely  that  these  conclusions  were  arrived  at  under  the 
idealized  conditions  of  the  simulations  scenarios.  The 
second  caution,  also  previously  stressed,  is  that  the  actual 
magnitudes  of  the  differences  between  position  estimates  has 
been  ignored.  It  is  conceivable  that  while  one  method 
always  produces  a  better  estimate  than  another,  the  differ¬ 
ence  between  any  two  position  estimates  is  acceptably  small. 

Lastly  there  is  the  caution  that  the  h.L.S.  model 
involves  a  complicated  iterative  procedure  whicn  uses 
considerable  computer  time.  It  is  probably  too  slow  a 
procedure  for  use  with  'real  time'  analysis  during  the 
execution  of  tracking  runs.  For  real  time  tracking  the 
NAV Y_A  method  currently  in  use  is  probably  preferrarle  cue 
to  its  simple  computations. 

However,  for  post  run  analysis,  and  also  possibly  for 
calibration  of  the  hydrophone  arrays,  the  h.L.S  method  is 
recommended  as  being  a  mere  exact  and  more  robi st  position 
estimator  than  those  methods  currently  in  use. 


I 
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C.  RECOMHENDATIONS  FOR  FOTOBE  INQUIRY 

Several  recommendations  have  already  been  made  for  work 
needed  to  estimate  the  effect  of  suitable  independent  vari¬ 
ables  on  both  timing  errors  and  the  bias  in  certain  methods. 
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In  addition  there  exist  at  least  two  other  areas  ior 
possibly  fruitful  investigation. 

The  first  area  concerns  the  interplay  netws en  methods. 
Specif ically,  the  binomial  comparisons  of  Chapter  V  show 
that  even  the  worst  methods  are  berter  than  each  of  the 
other  methods  at  least  part  of  the  time.  Hence  it  is 

possible  that  the  best  method  overall  would  be  a  suitable 
combination  of  methods,  wherein  each  method  is  used  where  it 
is  most  effective.  For  example,  even  thougn  the  M.L.3. 
method  has  been  indicated  as  the  best  method  for  any 
randomly  selected  position,  it  may  be  consistei  tly  outper¬ 
formed  by  another  method  under  certain  circumstances,  suer* 
as  extremely  high  or  low  elevation  angles. 

The  second  area  for  possible  work  addresses  the  question 
of  how  to  next  improve  upon  the  existing  model,  s.  It  is 
herein  suggested  that  the  next  improvement  in  modelling 
would  be  a  method  which  is  based  on  a  linear  velocity 
assumption.  As  figure  2.1  demonstrates,  a  linear  velocity 
profile  is  a  reasonable  approximation  for  most  depths.  This 
would  be  the  next  logical  step  above  the  constant  velocity 
assumption  which  is  associated  with  the  spherical  models. 
Most  of  the  mathematical  basis  for  such  a  model  in  contained 
in  Appendix  A.  Possibly  a  suitable  set  of  equations  could 
be  developed  involving  the  hydrophone  times  and  reflecting 
the  linear  velocity  assumption.  If  so,  the  least  squares  or 
maximum  likelihood  techniques  might  provide  useful  results. 
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APPENDIX  a 

LINE AS  VELOCITY  PROFILE  THEORY 

All  compu tations  in  this  study  are  tale  under  the 
assumption  that  the  sound  velocity  is  directly  related  to 
depth  in  a  linear  manner,  and  is  known  exactly.  Under  these 
circumstances  many  closed  form  results,  not  otherwise  avail¬ 
able,  can  te  obtained  and  used  in  those  computations. 

Suppose  that  the  velocity  profile  is  given  bf 

V  =  VO  +  V 1  *  z 

where  VO  and  VI  are  known  constants,  and  z  is  the  depth 
variable,  measured  down  from  the  water's  surface.  In  this 
case  it  is  known  [Bef.  4]  that  the  path  of  a  soi  nd  ray  from 
a  signal  source  to  a  hydrophone  will  have  the  shape  of  an 
arc  of  a  circle.  The  center  of  that  circle  will  be  some¬ 
where  above  the  surface  of  the  water.  The  vertical  place¬ 
ment  of  that  center  is  determined  by  tne  value  of  z  at  which 
the  speed  of  sound  equals  zero  (see  figure  A.1).  Although 
that  depth  is  negative  and  is  not  really  a  depth  at  all,  it 
nevertheless  has  geometric  meaning. 

Consider  the  vertical  plane  containing  the  circle 
center,  the  sound  source  and  the  hydrophone.  „et  h  be  the 
variable  which  measures  the  horizontal  position  in  that 
plane.  Let  (h ,  z  )  =  (a  1  ,a2)  be  the  position  of  the  hydrophone, 
and  let  (h,  z)  =  (p  1,  p2)  be  the  position  of  the  sound  source. 
Then  C2  given  by  (A.1)  is  the  Z  coordinate  of  the  circle 
center. 


What  must  be  found  is  Cl,  the  a  coordinate  o:  the  circle 
center,  and  r,  the  radius  of  the  circle.  To  solve  for  tuese 
values,  the  equation  of  the  circle  is  used,  evaluated  at  tne 
the  two  known  locations,  yielding  (A.  2)  . 
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Figure  A.1  Circular  Hay  Paths  for  a  Linear  Velocity  Profile 
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The  left  hand  sides  of  the  two  equations  of  (A. 2)  can  be 
equated,  and  solved  for  Cl,  leading  to  (A. 3). 
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(A. 3) 


When  this  value  of  Cl  is  substituted  inti  the  first 
equation  of  (A.  2),  then  the  result  is  given  by  (A.  4). 


The  circular  arclength  between  the  hydrophone  and  the 
sound  source  is  easily  computed.  Unfortunately  the  velocity 
of  sound  does  not  stay  constant  along  that  arc.  Therefore 
in  order  to  determine  the  amount  of  time  (T)  :  equired  for 

the  sound  ray  to  travel  the  ray  path,  the  effect  of  the 
velocity  must  be  integrated  along  the  arc.  This  is  dene  by 
(A. 5),  where  S  is  the  arclength  between  the  two  points,  and 
V  (s)  is  the  speed  of  sound  as  a  function  of  position  on  the 


Vis) 


(A. 5) 


In  (A. 5)  the  sound  source  position  corresponds  to  an 
arclength  of  s  =  0,  and  the  hydrophone  position  correponds 
to  arclength  s  =  S.  In  [Ref.  4]  this  integral  is  shown  to 
be  equivalent  to  (A.  6). 


d  a 

cos  i  a ) 


(A. 6) 


In  (A. 6)  AO  is  the  angle  of  elevation  of  the  ray  path  at 
the  sound  source,  and  A1  is  the  elevation  ai gle  at  the 
hydrophone.  The  antiderivative  in  (A.  7)  can  be  used  to 
solve  (A. 6)  ,  leading  to  the  ray  path  transit  time  expression 
in  (A .8)  . 


f- 


s  la) 


I  +  sin  la) 
cos  la) 


(A. 7) 


cos  ( A  x  '1  +  sin(A0)y 

AOV  V1  +  Sln  ( 1  )  / 


cos  ( 


(A. 8) 


If  the  elevation  angle  at  any  point  along  the  arc  is 
denoted  ty  A,  then  {A. 9)  relates  tne  angle  to  thj  derivative 
of  z  with  respect  to  h  along  the  ray  path. 


tan  ( A ) 


(A. 5) 


Implicit  differentiation  of  the  equation  of  the  circle 
yields  (A.  1 0)  as  another  expression  for  tae  same  derivative. 


h  -  c . 


(A. 10) 


Equating  these  two  derivatives  leads  to  (1.11)  which 
relates  the  elevation  anal e  and  the  position  (h,z)  on  the 


tan  ( A ) 


(A.  11) 


From  (A.  11)  a  simple  geometric  argument  produces  the 
equations  in  (A.  12)  . 


2  -  C. 


h  -  C. 
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(A.  12) 


First  let  A,  h  and  z  be  equal  to  (Al,a1,a2)  in  (A. 12), 
and  then  let  them  equal  {A0,p1,p2).  Then  substitute  both 
expressions  into  (A.  8).  The  result  is  (A.  13),  the  desired 
expression  for  the  ray  path  transit  time  in  te  rms  of  the 
positions  of  the  sound  source  and  the  hydrophone. 


T 


(A.  13) 


To  summarize,  if  a  ray  path  is  to  terminate  i  t  the  three 
dimensional  position  (X1,Y1,Z1),  and  the  sound  source  is  at 
(X2,Y2,Z2)/  then  perform  the  following  steps  in  order  to 
calculate  the  exact  ray  path  time  and  elevatioi  angle  that 
would  correspond  to  a  linear  velocity  profile: 
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APPENDIX  B 

PARTIAL  DERIVATIVE  FORMULAE  FOR  NEKTON' S  METHOD 


The  central  tool  used  b y  Newton's  Method  in  the  develop¬ 
ment  of  the  M.L.S.  model  in  Chapter  IV  is  the  matrix  GP.  It 
is  the  matrix  of  first  order  derivatives  of  the  error 
expressions  gl,  g2,  g3  and  g4.  Those  derivat.  ves  are  set 
forth  in  this  appendix. 

If  X  =  (01,02,03/0),  then  the  error  functions  are 

defined  by  (B.  1)  , 


g  .  (  x  ) 

i  — 


(x) 


T  „  -  D  M  /  V 

U  -  — - - 

1  -  N 


where  the  values  Ki,  M  and  N  are  the  f unctioi s  given  by 

(3.2)  . 
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Recall  that  GP  is  the  matrix  shown  in  (3.3) 


GP  (C  fC0,C3,U)  = 


Bgi 

Bgi 

dg. 

ag2 

a  *2 

a =2 

3  =  3 

Sg3 

dCl 

ac2 

SC3 

ac4 

Sq4 

ag4 

Bg4 

ag4 

Bc2 

dc3 

(3.3) 


Then  given  the  definitions  above,  the  following  deriva 
tives  are  the  result  of  straight  forward,  though  tedious 
differentiation,  and  are  offerred  without  detailed  proof. 
LEMMA  1 

V.  V  K ,  ( T ,  -  K ,  )  +  c  .  T  .  U  D 


LEMMA  2 


LEMMA  3 


LEMMA  4 


FORMULA  1 


The  first  three  diagonal  elements  of  GP  are  tne  deriva¬ 
tives  of  each  gi  with  respect  to  Ci  (i=1,2,3)  ai  a  are  given 
hy  (5.8). 


(3.8) 


FORMULA  2 

The  off  diagonal  elements  in  the  first  three  rows  and 
columns  of  GP  are  the  derivatives  of  each  gi  with  respect  to 
Cj,  and  are  given  by  (B .  9)  . 

i  -  1,2,3 
j  =  1,2,3 

j  *  i 


FORMULA  3 

The  derivative  of  each  gi  with  respect  to  U  is  given  by 

(B.  1  0) . 


T . M (VU  -  DC . )  +  VK . (T .  - 

l  l  li 


-  »  3  M 


(3.  10) 


V  M  “  K  .  /  K  . 

i  V  1 


FORMULA  4 

The  first  three  elements  of  the  last  row  of  G?  are  the 
derivatives  of  g4  with  respect  to  each  Ci,  and  ire  given  by 


(B  .  1  1 )  . 


3^4 


V  (1  -  N)  " 


3m 

&ci 


(l-W) 


(3.  11) 
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APPENDIX  c 

ON  IFOR  I!  SAMPLES  ON  A  TRONCAIED  HEMISPHERE 


The  single  array  simulation  of  Chapter  V  required  a 
random  sample  of  positions  in  space  uniformly  distributed 
over  the  surface  of  a  truncated  hemispnere.  Thi  hemisphere 
is  tc  be  of  radius  r  (3000  feet)  about  the  acoustic  center 
of  a  hydrophonic  array.  The  truncation  of  the  overall 
sphere  is  due  to  the  fact  that  the  upper  portion  of  the 
hemisphere  is  above  the  water  surface,  and  the  lower  half  is 
below  the  sea  bottom. 


w  =  r  cos (E)  dH 


Figure  C.1  Hemispherical  Geometry. 

Let  E  be  the  variable  denoting  angles  of  elevation  above 
the  horizontal.  Let  H  be  the  variable  denoting  horizontal 
azimuth  angles  about  the  center  of  the  hemisphere.  Consider 
a  small  piece  of  hemispherical  surface  area  bounded  by  the 


9  1 


elevation  angles  21  and  (ZI+dE),  and  by  the  azi  auth  angles 
HI  and  (HI+dK)  (see  figure  C.1).  If  W  is  the  horizontal 
width  of  that  piece,  then  5?  is  given  by  (C.1)  . 

W  =  r  dH  =  r  cos(E)  dH  (C.1) 

If  A 1  is  the  area  of  that  piece,  then  Al  is  i  pproximatea 
by  (C.2)  . 


A1  =  w  dE  =  r  cos  (E)  dH  dE  (C.2) 

Now  suppose  that  (  3  <  El  <  22  <  Eou  x  )  wnere 
Emax  is  the  elevation  angle  of  the  top  of  tne  truncated 
hemishpere.  Then  the  ratio  (A1/A2)  of  two  different  areas 
at  elevation  angles  El  and  E2  is  given  ny  (C.  3)  . 

A  r  cos (E^ )  dH  dE  cos(E,) 

A2  "  r  cos(E2)  dH  dE  =  "cosTeJT  ('-•3) 

If  nl  and  n2  are  to  be  the  (relative)  sample  sizes  from 
the  two  areas  Al  and  A2  respectively,  then  uniformity  of  the 
sample  requires  that  (C.4)  hold. 


AI  /  nl 


a2  /  n2 


(C.4) 


The  combination  of  (C.4)  and  (C.3)  implies  (1.5)  . 


cos  (E_, ) 


(C.5) 


letting  El  =  0,  then  the  relative  sample  size  n2  for 
area  A2  is  given  by  (C.  6)  ,  where  n  is  the  relative  sample 
size  at  the  base  of  the  hemisphere. 


n2  =  n  cos  (E0 ) 


(C  .  6 ) 


Now  the  differential  probability  of  drawii g  a  random 
position  that  has  elevation  angle  E  is  given  by  (C.7),  where 


N  is  the  total  sample  size  to  be  drawn  from  the  surface;  on 
the  sector  bounded  by  the  azimuth  angles  HI  and  (HI+dH). 


f„(E)  =  — 

L  N 


•cos(E) 


(C  •  7) 


Therefore,  if  K  is  defined  to  ne  the  constan:  ratio  n/N, 
then  the  corresponding  cumulative  distribution  is  given  by 
(C.  8)  . 


F^E) 


f  K  cos (e )  de 

J0 


K  sin (E) 


(C.3) 


For  the  highest  elevation  angle,  Emax,  th»  cumulative 
distribution  function  must  egual  one,  so  that  (C.9)  holds. 


f^e  ) 

E  max 


K  sin (E  ) 
max 


(C.9) 


As  a  result,  the  constant  K  is  determined  by  (C.10)  . 

v  _  l _ _ 

(C.10) 


sm(E  ) 
max 


Therefore  the  complete  cumulative  distribution  function 
is  given  by  (C .  1  1)  . 


F„  (E) 


sin (E) 


sxn(E  ) 
max 


(C.  11) 


Now  the  inverse  probability  transform  can  be  used.  If  U 
is  a  uniform  random  variable  on  the  interval  (0,1),  then  let 
E  be  given  by  (C.12). 

E  =  aresm  ( - ^ 

\ sm  (E  )J  (C.  1z) 

max 

Now  choose  an  azimuth  angle  H  randomely  ai  d  uniformly 
over  the  interval  (0,27T).  Then  (X,Y,Z)  given  by 

X  =  r  cos(fi)  cos  (E) 

Y  =  r  sin(H)  cos  (E) 

Z  =  r  sin(E) 


is  the  corresponding  position  in  spherical  roordinates 
That  position  is  a  random  position  drawn  from  a  popuiatio 
or  positions  uniformly  distributed  across  tne  si rface  of 
truncated  hemisphere  with  maximum  elevation  angle  Emax. 
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APPENDIX  D 

SINGLE  ARRAY  SIMULATION  COMPUTER  PEOGREB 


THIS  FORTRAN  PROGRAM  SELECTS  A  RANDOM  SAMPLE  OF 
SIZE  M  OF  3-DIMENSIONAL  POSITIONS  ALL  AT  A  RA”GE 
OF  F.  FEET  FROM  THE  ACOUSTIC  CENTER  OF  A  h'YDROPHONIC 
ARRAY.  THE  ARMS  OF  THE  ARRAY  ARE  ASSUMED  TO  EE 
ALIGNED  WITH  THE  COORDINATE  SYSTEM  OF  THE  RANGE. 

AND  IS  IN  A  POSITION  GIVEN  BY  THE  VECTOR  A. 

THE  SPEED  OF  SOUND  PROFILE  IS  ASSUMED  LINEAR. 


GIVEN  BY 


V  =  VV  (1) 


VV  (2)  *  DEPTH 


FOR  EACH  OF  THE  PCSIIONS  SELECTED.  EXACT  HYDROPHONE 
TIMES  ARE  CALCULATED  USING  THE  SUBROUTINE  TCOMP,  AND 
RANDOM  ERRORS  ARE  ADDED  TO  THE  EXACT  TIMES.  THE 
EP.ROP  DISTRIBUTION  IS  NORMAL,  WxTH  MEAN  ZERO  AND 
STANDARD  DEVIATION  SDEV. 

THE  TIME  VALUES  ARE  THEN  USED  TO  ESTIMATE  AN  APPARENT 
POSITION  BY  TWO  DIFFERENT  METHODS,  USING  APPROPRIATE 
SUBROUTINES.  THE  RESULTING  APPARENT  POSITIONS  ARE 
THEN  PROCESSED  BY  THE  SUBROUTINE  TRACE  TO  OBTAIN 
THE  CORRESPONDING  ACTUAL  POSITION  ESTIMATES.  THE 
RESULTING  POSITION  ESTIMATES  ARE  THEN  COMPARED  TO  THE 
ORIGINAL  TRUE  POSITION  TO  SEE  WHICH  METHOD  PRODUCED 
THE  MOST  ACCURATE  EST  IMA  IE  . 

THIS  IS  DONE  FOE  ALL  POSSIBLE  PAIRS  OF  THE 
SIX  METHODS. 

INTEGER  M,I, J, K.N. MET H. M ETH 1 , ME TH2, N YES . NY EST , NNOISE 


TEST  (10  00)  ,  ? T  ( 1  00  0 , 3)  .  DII 
DIFI  (1000[,AC(3[,TC  (  10 00 > 
R  ,  Z  ,  SDEV, PI, AZ  (10  00^,51^ 


DOUBLE  PRECISION  DIFT  M  0  00),  AC  (3)  .  ' 
DOUBLE  PRECISION  R ,  Z,  SDE  V ,  PI ,  AZ  (1  0 
DOUBLE  PRECISION  V,FK ,PKI, PHIMAX, S 

SET  INITIAL  VALUES  OF  SAMPLE  SIZE, 
AND  ERROR  STANDARD  DEVIATION 


000,4) 
SE. FEAC: 
F (1000) 

'10  00) 

t (loco,; 


a  =  iooo 

R  =  3000. DO 
SDEV  =  3. 4  64D-  5 
FORMAT  (2  X  ,  1  NOISE  STD. 


•  ,  FI  5.  10) 


FORMAT  (2  X .  '  NOISE  STD.  DEV.  =  ',F15. 

WRITE  (6,  1 1  SDEV 
WRITE  (6.99  9) 

SEED  =  93  1  947.  0D0 
PI  =  2. DO*  DARCOS(O.DO) 

SET  ARRAY  POSITION  AND  LINEAR  PROFILE 
FOR  SPEED  OF  SOUND. 


A  (1)  =  0.  DO 
A  2  =  0  .  DO 

A  ] 3)  =  1  3  00.  DO 
AC  (1)  =  A  7  11  + 

AC  (2)  =  A  (2)  + 

AC  (3)  =  A  (  3)  - 


15.0 
15.  0 
15.0 
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AD-A152  851  ALTERNATIVE  MODELS  FOR  CALCULATION  OF  ELEVATION  ANGLES 
AND  RAV  TRANSIT  TI.  .  OJj  NAVAL  POSTGRADUATE  SCHOOL 
NONTEREV  CA  C  D  MAIN  SEP  84 
UNCLASSIFIED  F/G  17/1  NL 


MICROCOPY  RtSOLUTION  USI  CHART 


non  nono  non  nnno  non 


10 


20 


100 


122 

11  1 


30  1 


100 


V  (1)  =  4840. 7D0 

V  (2)  =  33  14.  D-5 
=  VV  (1 )  *77  (2)  *A  (3) 

SET  UPPER  LIMIT  ON  ELEVATION  ANGLE,  IE  NEC2  SSAF.Y 


VV 

VV 

V 


I?  (  R  .GT.  A  ( 31  )  GO  TO 
PHIMAX  =  1.2566D0 
GO  TO  2  0 

PHIMAX  =  DAP.S1N  (A  (3)  /R) 
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GENERATE  RANDOM  VALUES  FOR  ANGLES  OF 
ELEVATION  AND  AZIMUTH. 


FK  =  1.D0/DSIN  (PHIMAX) 
CALL  GGUBS  (SEED,M,RAZ' 
CALL  GGUBS  (SEED,M,REL 


CONVERT  ANGLES  TO  3-D  COORDINATES 
DO 


111=  1,M 
EL  (I)  =  REL  (I) 

AZ  (I  =  RAZ  if 
PHx  =  D  ARSIN  (EL  (I 
P 
P 


11  CONTINUE 


FK) 

1.1)  =  R*DCOS  (PHI)  *DCCS  (AZ  (I)  *?I*2.D0  * 

1.2)  =  B*ECOS  (PHI)  *DSIN(AZ  (I)  *?I*2.D0) 

1.3)  =  A ( 3) -R*DS IN  (PHI) 


COMPUTE  EXACT  HYDROPHONE  TIMES  UNDER  THE  LC NEAR 
VELOCITY  PROFILE  ASSUMPTION. 


CALL  TCOMP  (P  ,  AC,  M,  VV,  TC) 
CALL  TCOMP  ( P ,  A, M ,  VV,  T) 


GENERATE  AND  ADD  ERRORS  TO  TIMES. 


NNOISE  =  4  *M 

CALL  GGNML  (SEED, NNOIS  E, STUFF) 

DO~1 1  II  =  1 , M 
DO  122  J  =  1,4 

STUFF (K) 

NCISE*S  DEV 
=  T  (I,  J)  ♦NOISE 


NOISE  = 
NOISE  = 


OTi 


CONTINUE 

CONTINUE 


FORM  AT  (2  X  ,  » 
WRITE  (6,3  0  1) 
(6, 


WRITE 


(  SAMPLE  SIZE  =  ’ , I j 


9  9  9) 


TITO  CUTER  LOOPS  RUN  THROUGH  ALL  PAIRS  OF  HE  SIX 
PCSSI ELS  POSITION  ESTIMATING  METHODS. 


DO  1111  METH1  =1,5 
KMETH  =  METH1  ♦  1 
DC  1222  METH2  =  KMETH, 6 

IF  (METH1  .EC-  METH2)  GO  TO  1222 


FORMAT  (2X.  ’METHODS  USED  ARE 
WRITE  (6,  200) 

METH  =  METH1 


’  ) 


nnno  noon  nfinnn  nnn 


CALL  SUBROUTINES  TO  IMPLEMENT  THE  METHODS 
IF  (METH  .HE.  1)  GO  TO  131 


(METH  .  NE.  1)  GO  10  131 
CALL  POSNAV  (TN,V,  M, PEST, TEST) 
FORMAT  (2X, '  NAVY, 

WRITE  (6,20  1 ) 

GO  TO  150 

(METH  . NE.  2)  GO  1 0  132 
CALL  POSNVC  (TN  ,7,  M, PEST, TEST) 
FORMAT  (2X,'  NAVY, 


UNCORE E:  TED’) 


FORMAT (2X,'  '  '  '  *  '  NAVY,  CORRECTED') 

WRITE  (6  .202) 

GO  TO  1 50 

IF  (METH  .NE.  3)  GO  TO  133 

CALL  POSLS  (TN,  V,M  ,  PEST,  TEST) 

FORMAT  (1  8X  ,  'LEAST  SQUARES,  JNCOR  RECTED  ' ) 
WRITE (6,203) 

GO  TO  150 

IF  (METH  .  NE.  4)  GO  TO  134 

CALL  POSLSC  (TN  ,V,  M, PEST, TEST) 

FORMAT (18X, 'LEAST  SQUARES,  CORRECTED') 

WRITE (6,204) 

GO  TO  1 50 

IF  (METH  .  NE.  5)  GO  TO  135 

CALL  POSMLP  (TN,V.  M,  PEST,  TEST,  VAR) 

FORMAT ( 18X,  'MAX.  LIKELIHOOD,  PLANAR’) 

WRITE  (6,205) 

GO  TO  150 

IF  (METH  .  NE.  6)  GO  TO  136 

CALL  POSMLS  (TN , V, M, P EST.TES T, V AR) 

FORMAT (18x,  'MAX.  LIKELIHOOD,  SPHERICAL'  * 
WRITE  (6,206) 

GO  TO  1  50 

WRITE  (6,13  7)  METH1.METH2 

FGRMATj2X, 'CANT  FIND  METHODS  ',14,'  AND/OR  ’,14) 

GO  TO  fOOfi 

THE  APPARENT  POSITION  ESTIMATES  ARE  FTLAIIVE  TO 

THE  ARRAY  CENTER,  AND  MUST  EE  TRANSLATED  TO 

THE  TRACKING  RANGE  COORDINATE  SYSTEM. 


150 

DO  144  I  = 

1,M 

PEST  (I, 

1  = 

PEST 

(I 

PEST  I, 

2  = 

PEST 

(I 

144 

PEST  I , 
CONTINUE 

3  = 

A  (3) 

CONVERT  APPARENT  POSITIONS  TO  ACTUAL  ESTIMATES  BY 
RAY  TRACING  PROCEDURE. 

CALL  TRACE  (P SST, PT,  T EST , A , M, VV) 

CALCULATE  DIFFERENCE  BETWEEN  THE  1ST  POSITION 
ESTIMATE  AND  THE  TRUE  POSITION. 


(  METH  .NE.  MET  HI  )  GO  TO  160 
DO  155  I  =  1 , M 


DIF  (I)  =  0. DO 

DIFT(x)  =  DABS  (TEST  (I) -TC  (1,4)  ) 

DO  156  J  =  1,3 

DIF  (  I)  =  DIF  (  I)  +  (PT  (I,J) -P  (I,J)  )  **l 
CONTINUE 


CONTIN 

CONTINUE 


nnnnn 


ME  1H  =  MET  H  2 
GO  TO  1  30 

CALCULATE  DIFFERENCE  BETWEEN  THE  2ND  POSITION 
ESTIMATE  AND  THE  TRUE  POSITION,  AND  COMPARE  IT 
TO  THE  1ST  POSITION  DIFFERENCE. 


NYES  =  0 
NYEST  =  0 
DO  166  I  =  1 , K 
DO  167  J  =  1.3 

DIF  (I)  =  DIF  (I)  -  (PT(I,J)-P(I,J))  **2 

CONTINUE 

D I  FT  ( I)  =  DIFT  (I)  -  D  ABS  (TEST  (I) -TC  (1,4)  ) 
IF  (  DIFT(I)  .GT.  0.  DO  )  GO  TO  165 
NYEST  =  NYEST  +  1 
IF  DIF(I)  .GT.  0  .  DO  )  30  TO  166 

NYES  =  NYES  +  1 
CONTINUE 


FRACT  =  (DFLOAT  (NYES)  /DFLOAT(M)  ) 
FORMAT  (2X,*  FRACTION  OF  TIME  M^Tr 


THOD  #1  IS  :  LOSER* ) 


,  F7. 31 
, FT.  3] 


C 

1222 

11  1  1 

C 

993 

9S9 

1000 


FORMAT  (2  X  ,  *  IN  POSITION  :  *,F7.3) 

FORM  AT  (2  X  .  '  IN  TIME  .  \F7.3) 

WRITE  (6,  3  0  0) 

WRITE  (6,303)  FRACT 

FRACT  =  (DFLOAT  (NYEST)  /D FLOAT  (H )  ) 

WRITE  (6,  jO  2)  FRACT 
WRITE  (6,  999 
WRITE  6,  998) 

WRITE  (6,  999 

CONTINUE 

CONTINUE 

FORMAT  (IX,  *=  =  =  =  ==  =  =  =  =  =  =  =  ==  =  =  =  •) 
FORMAT  (2  X  ,  •  •) 

STOP 
END 
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APPENDIX  E 

DOUBLE  ARRAY  SIMULATION  COMPUTER  PROGRAM 


THIS  FORTRAN  PROGRAM  SELECTS  A  RANDOM  SAMP',  E  CF 
SIZE  M  OF  3-DIMENSIONAL  POSITIONS  IN  A  BOX  RUNNING 
CROSSWAYS  BETWEEN  TWO  HYDROPHONIC  ARRAYS.  THE  30X 
HAS  DIMENSIONS  GIVEN  BY  THE  VECTOR  BOX.  THE  ARRAYS 
ARE  SEPERATED  EY  7500  FEET.  AND  HAVE  COORDINATES 
GIVEN  BY  THE  VECTORS  A1  AND  A2.  THE  ARMS  OF  THE 
ARRAYS  ARE  ALIGNED  WITH  THE  RANGE  COORDINATE  SYSTEM 


THE 

AND 


SPEED  OF 
IS  GIVEN 


SOUND 

BY 


PROFILE 
V 


E  IS  ASSUMED  TO  3E  LINEAR, 
=  VV(1)  +  VV  (2)  *  DEPTH  . 


FOR  EACH  OF  THE  POSITIONS  SELECTED,  EXACT  TIME  VALUES 
FOR  SOUND  WAVE  ARRIVAL  AT  THE  HYDROPHONES  OF  EACH 
ARRAY  ARE  COMPUTED  BY  THE  SUBROUTINE  TCOMP.  THEN 
RANDOM  ERRORS  ARE  ADDED  TO  ALL  TIMES.  THE  ERROR 
DISTRIBUTION  IS  NORMAL,  WITH  MEAN  ZERO  AND  STANDARD 
DEVIATION  SDE  V . 


THE  TIME  V 
POSITIONS 
A  PFROPRI AT 
DIFFERENT 
POSITIONS 
ESTIMATES 
CF  THE  DIF 
THE  TWO  DI 
METHOD  PRO 


ALUSS  ARE  THEN  USED  TO  ESTIMATE  IPPAEENT 
BY  TWO  DIFFERENT  METHODS,  USING  THE 
E  SUBROUTINES.  EACH  METHOD  PRODUCES  TWO 
ESTIMATES,  ONE  FOR  EACH  ARRAY.  THE  APPARENT 
ARE  THEN  TR ANSLATED  TO  ACTUAL  POSITION 
BY  THE  SUBROUTINE  TRACE.  THEN  I  HE  LENGTH 
FERENCE  VECTOR  FOR  EACH  METHOD  IS  COMPUTED. 
FFERENCE  LENGTHS  ARE  COMPARED  TO  SEE  WHICH 
DUCES  POSITION  PAIRS  IN  CLOSEST  AGREEMENT. 


THIS  COMPARISON  IS 
THE  SIX  METHODS. 


DONE  FOR  ALL  POSSIBLE  Pi  IRS  OF 


INTEGER  M , I, J, K,METH,  METHl , METH2, NYZS,NNOI3  E ,  N  AR R  AY 
DIMENSION  RP  (  1  000)  ,STUFF  (4 0 00)  ,  NAM  1  (6)  ,  NAM 2  (6) 

1000,4) 


DOUBLE  PRECISION  A  1  (3  )  ,  A  2  { 3)  ,  T1  (1  0  00 , 4)  ,  T2( 

DOUBLE  PRECISION  V AP (  lOOOi , NOIS S, FT - - 

DOUBLE  PRECISION  TEST 

D0U3LE  PRECISION  SDEV 

DOUBLE  PRECISION  T(1, 


IP.  (  100  0)  .NOISE.  FR  ACT,  XYZ.  DEPTH 
;ST  (1000)  ,?  1  (  1000,  3)  ,  ?2(f  0  00,4) 
)EV,SEED,30X(3)  ,p5Si(10J0,3)  ,A  (3) 
[10  00,4)  ,V,VV  (2)  ,  DIF  ( 1 00  0) 


DATA  NAM1 

(1 

,  N  AMI 

[2 

,  NAM1 

l3) 

DATA  NAM  1 

,4, 

,  NAM  1 

5 

,  NAM  1 

6 

DATA  NAM 2 

1 

,  N  AM2 

’2 

,  NAM 2 

3; 

DATA  NAM  2 

4 

,  NAM2 

(5 

,  NAM  2 

[o) 

/’ NAVY’ , *  NAVY* 


/'L 

/» 

/’COi 


t  n  A  »  A  f 

.  '  ,  *M.  L.  •  , 

' ,'CORR'  * 
R'  » 


1  PL  NR 


’L.S. ’/ 
•  M.  L.  V 


.  V 

'  SPHR  ’  / 


INITIALIZE  VALUES  FOE  SAMPLE  SIZE, 
AND  ERROR  STANDARD  DEVIATION. 


RANGE,  DEPTH 


M  =  1000 
SDEV  =  3.464D-5 
DEPTH  =  1300. DO 
FORMAT  (2  X,  '  NOISE 
WRITE  (6,  if  SDEV 
WRITE  (6.  999) 

SEED  =  9 34  7.  6  D  0 


STD.  DEV. 


=  i 


F 1  5.  10) 


99 


-1 


SET  0?  ARRAY  POSITIONS,  AND  SOUND  VELOCITY  PROFILE 
[1)  =  -3750  . DO 


(ij  =  0.D0 


=  DEPTH 
=  3750. DO 
=  0.D0 
=  DEPTH 
=  4640. 7D0 
=  33  14.D-5 
=  V V  (1 )  +  VV  (2)  *DEPTH 


FORMAT  (21.  • 
WRITE  (6,30)  M 


{  SAMPLE  SIZE  = 


WRITE 

DRAW 


ii'99  9) 

UNDERWATER  BOX  FOR  RANDOM  POSITIONS 


)  ’) 


BOX 

BOX 

BOX 


il}l 


1 000. DO 
5000.  DO 
DEPTH 


GENERATE  RANDOM  POSITIONS  IN  BOX 


DO  1  1  J  =  1,2 

CALL  GGUBS  (SEED,  M  , RP ) 

DO  22  I  =  1,M 
XYZ  =  RP  (I) 

P1(I,J)  =  (XYZ-  0.  5D0)  *BOX  (J) 
CONTINUE 
CONTINUE 

CALL  GGUBS  (SEED, M,RP) 

=  i,M 


DO  33  I 
XYZ  = 


CONTI 


PI  (1,3?  VI  YZ*BOX  (  3) 


COMPUTE  EXACT  HYDROPHONE  ARRIVAL  TIMES,  AND  ADD 
RANDOM  ERRORS  TO  TIME  VALUES. 


NNOISE  =  M  *4 

CALL  TCOMP  (P1,A1,M,VV  ,T  1) 

CALL  GGNML  (s£eD,NNUI  SE,  STUFF) 


K  =  1 
DO  44  I  =  1 , M 

DO  45  J  =  1 


noise  =  Stuff (K ) 

NOISE  =  NOISE* S  D 


T1  (I  -  J) 
K  =  1 

CONTINUE 
CONTINUE 


£*  V 

=  T 1  (I,  J)  ♦ 


NOISE 


CALL  TCOMP  (Pi  ,A2,  M,VV  ,T2) 

CALL  GGNML  ( S EEE, NNO I SE, ST UFF) 

DO  55  I  =  1 , M 
DO  56  J  =  1,4 

NDISE  =  STUFF  (K) 

NOISE  =  NCIS  E*S  DEV 
T2  (I,J)  =  T2  (I,  J)  *  NOISE 


nnnoo  n  no  nnn 


SET  Ur  LOOPS  TC  HUN  THROUGH  ALL  PATHS  Or  M2THODS 


DO  1111  METE  1  =  1,5 
KMETH  =  MET  HI  ♦  1 


DO  1222  METH2  =  KMETH .6 

IF  (METH1  .EC.  MET  n2)  GO  TO  1222 


FORMAT  (2X,  'METHODS  COMPARED  ARE  1. 
FOEMAT  (2  X  ,  *  AND  2. 

WRITE  (6/  10)  NAM1  (METH  1)  ,  NAM2  (METH  1) 
WRITE  (6,  20  j  NAM1(METH2)  ,  NAM2  (METH2) 
WRITE  6,999) 


1.  •  ,  A  4,  2  X  ,  A4) 

2.  *,A4,2X,A4) 


METH  =  METH1 

NARRAY  =  1 

SET  UP  ARRAY  EEIfJG  USED 

i\l\  :il|3 

sJ3Uvi‘n« 

DO  67  J  =  1,4 

T  (I .  J)  =  T1  (I,  J) 

CONTINUE 

CONTINUE 

CALL  SUBROUTINES  TO  PERFORM  ESTIMATION  METiODS. 

IF  {METH  .  NE.  1)  GO  TO  131 

CALL  POSNAV  (T  , V,  M  , PEST,  TESI ) 

GO  TO  150 

IF  (METH  .NE.  2)  GO  TO  132 

CALL  POSNVC  (T,7,M  , PEST, TEST) 

GO  TO  150 

IF  (METH  .NE.  3)  GO  TO  133 

CALL  POSLS  (T,V,M,  PEST,  TEST) 

GO  TD  150 

IF  (METH  .NE.  4)  GO  TO  134 

CALL  POSLSC  (T,V,M  , PEST, TEST) 

GO  TO  150 

IF  (METH  .NE.  5)  GO  TO  135 

CALI  POSMLP  (T,V,M,  PEST,  TESI,  VAR) 

GO  TO  150 

IF  (METH  .  NE.  6)  GO  TO  136 

CALL  POSMLS  (T  ,  V,  M  ,?E ST,  TEST  ,  VA R) 

GO  TO  150 

WRITE  (6,  13  7)  METH1.METH2 

FORM  AT(2X, ’CANT  FIND  METHODS  ',14,'  AND/OR  ’,14) 

GO  TO  1000 


(  NARRAY 


2  )  GO  TO  152 


APPARENT  POSITION  ESTIMATES  ARE  IN  LOCAL  A?  RAY 
COORDINATES,  AND  MUSI  3E  TRANSLATED  TO  TRACKING 
RANGE  SYSTEM  COORDINATES. 


DO  144  I  =  1 , M 

PEST  (I,  1)  =  PEST  (I  ,  1) 
PEST  1,2)  =  PEST  (I  ,2  [ 

PEST  1,3)  =  A1  (3)  -  P 

CONTINUE 


♦  ai  (i; 

♦  A 1  2 
EST  (1,3 


CORRECT  APPARENT  POSITIONS  3Y  RAY  TRACING. 
CALL  TRACE  (P  EST,  P  1,  T  EST  ,  A  ,  M,  VV) 


noon 


GO  CN  TO  SECOND  ARRAY 


NARRAY  =  2 

:iiia 

A]3)  =  A  2  3 
DO  77  I  =  1,  M 


DO  78  J  =  1,4 

T(I,  J)  =  12  (I,  J) 

CONTINUE 
CONTINUE 
GO  TO  130 

TRANSLATE  2ND  ARRAY  APPARENT  POSITIONS  TO  i  ANGE 
COORDINATE  SYSTEM,  AND  THEN  RAY  TRACE. 


DO  145  I  =  1 ,  M 

PE  ST  (I,  1)  =  PEST  (I  ,  1)  +  A2  (1 ) 

PEST  ’1,2)  =  EEST(I,2)  +  A2  (2 
PEST  ’I,  3  =  A2(3)  -  PEST  (1,3) 

CONTINUE 


CALL  TRACE  (PEST, P2, T EST,A, M, VV) 

COMPUTE  DIFFERENCE  VECTORS  FOR  1SI  METHOD 

IF  (  MET  H  .HE.  MET  HI  )  GO  TO  160 
DO  155  I  =  1  ,M 
DIF (I )  =  0.D0 
DO  156  J  =  1,3 

DIF(I)  =  DIF  (I)  +  (PI  (I,  J) -P2  (I,  J)  )  **2 
CONTINUE 
CONTINUE 

GO  CN  TO  SECOND  METHOD 

METK  =  METH2 
NARRAY  =  1 
GO  TO  125 

COMPUTE  DIFFERENCE  VECTORS  FOR  2ND  METHOD  1  N D 
COMPARE  TO  DIFFERENCES  FOR  1ST  METHOD. 

NYES  =  0 
DC  166  I  =  1  ,  M 
DO  16  7  J  =  1,3 

DIF(I)  =  DIF  (I)  -  (P  1  (I,J) -?2  (I,  J)  )  **2 
CONTINUE 

IF  {  DI F  ( I)  .GT.  0  .  DO  )  GO  TO  166 
NYES  =  NYES  ♦  1 
CONTINUE 

WRITE  (6,  999) 

FRACT  =  (DFLOA  T  (NYES)  /DFLGAT(M)) 

FORMAT  (2  X,’  FRACTION  OF  TIME  METHOD  #1  ’) 


:  (2X,  ’ 
FORM  AT  (2  X  ,  ' 
WRITE  (6,  300) 
WRITE  (6,  3  10) 
WRITE  (6,  999) 
WRITE  6, 993) 
WRITE  (6,  99  9) 

CONTINUE 

CONTINUE 


1222 

1111 


IS  BETTER  IS 


•  ,F7. 3  -- 


FRACT 


FORMAT  (1  X, 
FORMAT  (2  X,  ' 
STO 


no  non  onn  nnnnnnnnnnnn  nn 


APPENDIX  F 

CO HPOTEE  SUBROUTINES  FOR  TIME  CALCULATION  ASD  RAYTRACIN G 


SSSSSSSSSSSS  SSSSSSSSSSSSSSSSSSSSSSSSSSSSS5SSSS5 SSSSSSSS 


SUBROUTINE  T CO  ME  ( P, A  ,  M,  7V , T) 

THIS  FORTRAN  SUEROUIINE  COMPUTES  THE  EXACT  TIME 
REQUIRED  FOR  A  SOUND  WAVE  TO  TRAVEL  FROM  ITS 
SOURCE  (VECTOR  P)  TO  THE  FOUR  H I  DROP  HONES  ON  AN 
ARRAY  WHOSE  ACOUSTIC  CENTER  IS  SPECIFIED  BY 
VECTOR  A. 

THE  METHODOLOGY  USED  IS  SET  FORTH  IN  APPENDIX  A 
OF  THE  THESIS.  THE  BASIC  ASSUMPTION  IS  ONE  OF 
A  LINEAR  VELOCITY  PROFILE,  WHOSE  COEFFICIENTS  ARE 
GIVEN  BY  THE  VECTOR  VV. 


22 

11 


44 

33 


INTEGER  M,I,J 
DOUBLE  PRECISION 


P  (1 
AA  ( 


DOUBLE  PRECISION 
SET  UP  COORDINATES  OF 
DO 


000 

4 


FOUR  HYDROPHONES  ON  THE  ARRAY 


111=  1,4 
DO  22  J  =  1,3 

AA  (I,J)  =  -15.  D 

J»TT  MfT 


CONTI 


AA 

AA 


CONTINUE 

AA  (1,3)  =  15. DO  * 
TINUE 


0  ♦  A  (J) 
A  (3) 


1,  1 

2,2 


=  15. DO  ♦  A 
=  15. DO  ♦  A 


HI 


A  A  (3 , 3;  =  -15. DO  *  A'( 
CALCULATE  QUANTITIES 
C2  =  -1.  DO*  (VV  (1)/VV( 

M  SO 


LOOP  THROUGH  THE 
DO  33  I  =  1 , M 
P2  =  P  (1,3) 

LOOP  THROUGH  TH 
DO  44  J  =  1.4 

PI  =  DS  Q&T  (  ( P  (I 
♦  ( 

Cl  =  P1**2+P2** 

♦  2. DO 

Cl  =  C 1 /  (2. D0*P 
R  =  DSQRT  (C1**2 
T  (I  ,  J)  =  ((AA  (J 


Cl,  C2 ,  R,  AND  THE  TI!  ES  T{*,4) 

2)) 

URCE  LOCATIONS 


FOUR  HYDROPHONES  FOR  E\  CH  SOURCE 


.1) -AA  (J.  1)  )  **2 
P  (1 ,2)  -A  A ( J, 2) )  **2) 
2-AS(J,  3)**$ 

*C2* (AA (J, 3) —  P  2 ) 

1). 


:’j 


/  (P2-C2) 


T(I,  J) 
CONTINUE 
CONTINUE 


*'((fi-Cl+P^i/(R-CI  )  ) 
=  (DLOG  ( T  (I  ,  J)  )  )  /V  V  ( 2) 


FETURN 


nnn  nno  n  nono  o  non  o  onnnnnononoo  on 


SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 
SUBROUTINE  TRACE  (P,  PT  ,  T  ,  A,  M,  V  V) 

THIS  FORTRAN  SUBROUTINE  TAKES  AN  APPARENT  POSITION 
ESTIMATE  ?  AND  CONVERTS  IT  TO  AN  ACTUAL  POSITION 
ESTIMATE  PT  3Y  RAY  TRACING.  THE  RAY  TRACING  ASSUMES 
A  LINEAR  VELOCITY  PROFILE.  WITH  COEFFICIENTS  IN  THE 
VICTOR  VV.  OTHER  INPUTS  ARE  THE  ARRAY  LOCATION 
VECTOR  A  AND  THE  RAY  TRANSIT  TIME  TO  THE  ACOUSTIC 
CENTER.  T.  THE  RAY  TRACING  LAYERS  ARE  DEL  FEET  THIC 
WITH  THE  FIRST  (BOTTOM)  LAYER  BEING  THE  Eli  ST  DEL 
FEET  IMMEDIATELY  ABOVE  THE  ARRAY. 


INTEGER  M,  I 

DOUELE  PRECISION  T  M 0  00)  ,?(1000,3)  ,PT(10QQ,  3),A(3) 
DGU3LE  PRECISION  DEL . H ,S A, C A, DEPTH , V , V 1 , VV (2) 
DOUBLE  PRECISION  DS,DT,DZ,DH 

DEL  =  25. DO 

LOOP  THROUGH  THE  M  APPARENT  POSITIONS 
DO  11  I  =  1,  M 

INITIALIZE  INCREMENTAL  VALUES  FOR  EACH  POSITION 

!r.Ds?sii!  is/w 1 ,21' 1  **2) 

D  T  =  0 .  DO 
DZ  =  0. DO 
DH  =  0. DO 

DEPTH  =  A  (3)  -  DEL/2. DO 
V  =  VV  (1)  +VV  12)  *DEPTH 


(2)  *DEPTH 


INNER  LOOP  :  PROCESS  DATA  UPWARDS,  LAYER  BY  LAYER, 
UNTIL  P.  AY  TRANSIT  TIME  IS  EXHAUSTED 

SA  =  DSQRT  (  l.DQ-CA**2) 

DZ  =  DZ  ♦  DEI 
DS  =  DEL/SA 
DT  =  DT  ♦  DS/V 
DH  =  DH  ♦  D S*CA 

IF  (  DT  -GT.  T  (I)  1  GO  TO  20 

DEPTH  =  DEPTH  -  DEL 
VI  =  VV  (1)  ♦  VV  (2)  *DEPTH 

CA  =  CA*(V1/V) 

V  =  VI 
GO  TO  10 

ADJUST  FOR  OVERSHOOTING  IN  LA  ST  LAYER 

DS  =  V*  (DT-T  (I)  ) 

DH  =  DH  -  C A*DS 

ADJUST  APPARENT  POSITION  TO  GET  ACTUAL  POSCTION 

?  T  (1,1)  =  fP  Cl,  11  -  A  f  11  3  MDH/H)  +  AM) 

PT  (I,  2)  =  P  1,2  -  A]2i  *  DH/H  *A  2 
PT  (I,  3)  =  A  (3)  -  ( DZ-5A*DS) 

CONTINUE 

FETURN 


nnnnnnnnnnnn  n  on 


APPENDIX  G 

COMPUTES  SUBROUTINES  FOR  POSITION  ESTIMATION  METHODS 


SSSSSSSSSSSSSSSSSSSSSSSSS  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 
SUBROUTINE  P0SNA7  (T,V,M,?,NT) 

INTEGER  M,  I 

THIS  FORTRAN  SUBROUTINE  IMPLEMENTS  THE  ORIGINAL 
UNADJUSTED  NAVY  APPARENT  POSITION  ESTIMATION  METHOD- 

INPUT  ARE  THE  HYDROPHONE  TIMES  T  FOR  M  SOUi  D  SOURCE 
POSITIONS,  AND  VELOCITY  V  AT  THE  ARRAY. 

OUTPUT  ARE  M  APPARENT  POSITION  ESTIMATES  ?  ALONG  WITH 
THE  CORRESPONDING  M  RAY  TRANSIT  TIMES  NT.  ALL  THE 
POSITIONS  ARE  REFERENCED  TO  THE  ACOUSTIC  C2NTER,  WITH 
THE  2  COMPONENT  BEING  MEASURED  UPWARDS  FPOM  THE  ARRAY. 

DOUBLE  PRECISION  T  ( 1 0  00 , 41  ,  P  (  1 0  00 , 3 )  ,  NT  ( 1  0)  0) 

DOUBLE  PRECISION  V,TC ,D, DISC, NU MEE, DENOM 

D  =  30. DO 
DC  11  I  =  1.  M 
NUMER  =  O.DO 
DENOM  =  O.DO 
TC  =  T  (I,  4 ) 

DO  22  J  =  1,3 
P(I,J)  =  V*V* 

NUMER  =  NUMER 
DENOM  =  DENOM 
22  CONTINUE 

NT  (I)  =  TC*DSQRT  (N  UMER/DENOM ) 

11  CCNTI NUc 
FETURN 


(TC*TC-T  (I.  J)  **2)  /{D*2.  DO) 
*  Pjfl,  J)**2 
+  (?(I,J)  ♦  1 5.  DO)  **2 


nnnnnnnnnnnnn  nn 


SSSSSSSSSSSSSSSSSSSSSSSSS SSSSSSSSSSSSSSSSSSSSSJ  SSSSS  SS3 SSS 


SUBROUTINE  P  OSNVC  (7,  V,  K#P,HT) 

THIS  FORTRAN  SOBF.OUTINE  IMPLEMENTS  THE  ADJUSTED  NAVY 
METHOD  N  AV  Y_  A  FOR  ESTIMATING  APPARENT  POSITION'S. 

INPUT  ARE  THE  HYDROPHONES  TIMES  T  FOR  M  SOJ ND  SOURCE 
LOCATION S ,  AND  THE  VELOCITY  V  AT  THE  ARRAY. 

OUTPUT  ARE  THE  M  APPARENT  POSITION  ESTIMATES  ?,  AND 
THE  CORRES  PONDING  M  RAY  TRANSIT  TIMES  NT.  ALL 
POSITIONS  ARE  P.EFSFENCED  TO  THE  ACOUSTIC  TENTER  OF 
THE  ARRAY.  WITH  THE  Z  COMPONENT  MEASURE  UPWARDS 
FROM  THE  ARRAY. 


INTEGER  M , I 
DOUBLE  PRECISION 
DOUBLE  PRECISION 

D  =  30. DO 
DO  1  1  I  =  1, M 
DCC  =  0 . DO 
TO  =  T  (1,4) 

DO  22  J  =  1  , 3 


T(1000,4),P(1000,3),NT(1}00) 
V,  D  ,  NUMEE  ,  DeN  CM ,  TO  ,  DCC 


P(I,J)  =  15. DO  ♦  (V*V*(TC*IC-T  (I,  J)  *«  2)  /  (D*2.  DO 
DCC  =  DCC  ♦  P  (I  ,  J)  **  2 

TTNfT  V. 


CONTINUE 
NUMER  =  0.D0 
DENOM  =  0.D0 
DO  44  J  =  1,3 


/DS  2RT  (DCC) 


CC  NT IN 
FETURN 
E 


CONTINUE 

Tllful  =  rC*£S£fiT(NUJ1Ei’'/D£N0!^ 


SSSSSSSSSSSS 


SSSSSSSSSS SSSSSSSSSSSS  SSSSSSSSSSSS SSSSSSi 

S'JE ROUTINE  POSLS  (7,  V  ,  M,  ?,  LSI) 

THIS  FORTRAN  SUBROUTINE  IMPLEMENTS  THE  LEAST  SQUARES 
PLANAR  METHOD  I.S. 


INPUT  ARE  THE  HYDROPHONE  TIMES  T  FOP.  M  SOUND  SOURCE 
POSITIONS,  AND  THE  VELOCITY  V  AT  THE  ARRAY. 

INPUTS  AND  OUTPUTS  ARE  THE  SAME  AS  FOE  THE 
SUBROUTINE  POSNAV. 


INTEGER  M.I.J 

DOUBLE  PRECISION  T  ( 1  0  00 , 4)  ,  P  (  1 0  00 , 3)  ,  LS  T  ( 1 )  00 ) 
DCUELE  PRECISION  V,DISC,TC 

DC  11  I  =  1,  M 

DISC  =  0 . DO 
DO  22  J  =  1,3 

PJI,J)  =  T  (I,  4)  ~T  (I,  J) 

DISC  =  DISC  +  P (I, J) **2 
CONTINUE 

LSTfl)  =  (T  (1,1)  +1  ( I, 2) +T  (I  ,3)  -  T  (1,4)  )/2.  DO 
DO  33  J  =  1,3 

P(I,J)  =  (V*LSI  ( I)  *P  (I ,  J)  )  /DSQP.T  (DISC) 
CONTINUE 
CONTINUE 


SSSSSSSSS2SS  333S5SS2SS33S  SSSSSSSSSSSooSSSSSSSSJ  SS3S3SS335: 
SUBROUTINE  PCSLSC  (T ,  V  ,  H  ,  P,  L3  T) 

THIS  FORTFAN  SUBROUTINE  IMPLEMENTS  THE  3IAS  ADJUSTED 
LEAST  SCUARES  PLANAR  METHOD  L.S.C.  FOR  ESTIMATING 
APPARENT  POSITIONS. 

INPUTS  AND  OUTPUTS  ARE  THE  SAME  AS  FOR  THE 
SUBROUTINE  POSTS. 


INTEGER  I  .  M.J.K.L 


1  H  iLVJ  X  *  l  1  m  m  m  ±J 

DOUBLE  ?  F.SCISIC N  T  { 1  0  00. 4)  ,P  {  1 0  00,  3)  ,  LSI  ( 1 )  00) 
DOUBLE  PRECISION  V  ,  A  (  4, 3 )  ,  D  (4  )  ,  DI  SC,  TD,  h  IS  0  R  ,  R 


:  R  *  I  ( 


SET  UP  HYDRO  PH  CNE  POSITIONS 
DO  44  J  =1,3 

DO  55  I  =  1 ,4 

A(I.J)  =  -15. DO 
CC  Nil  HUE 
A  (J. J)  =  1  5 .DO 
CC NT  I  HUE 

CAICULATE  ORIGINAL  L. S.  SOLUTION 
DO  11  I  =  1,  M 

DISC  =  0 . DO 
DG  22  J  =  1,3 


22 

CONTINU 
LSI  m 

DO  j3  J 

ccntInTj 

33 

R1SQF.  = 

K 

89 

cgnt!n!j 

TD  =  0. DO 
DO  77  J  =  1  ,  3 
D  fJ)  =  0.  DO 
DO  8  8  K  =  1,3 

D{ J)  =  D(J)  ♦  (?<I,K)  -  A  (J,K))  **2 

CONTINUE 

D  ( J )  =  D S CRT  (D  ( J }  ) 

DISC  =  DISC  +  (D  (4)  -D  (J)  )  **2 
TD  =  TD  ♦  D ( J) 

CONTINUE 


CALCULATE  BIAS  VECTOR  E 
DISC  =  D SQ RT  (DISC) 

DO  66  J  =  1,3 

E(J)  =  {  (TD-D{4)  )  *  (D  (4)  -  D  ( J)  ) /  (DISC* 2.  ,D0)  )  - 
CON  IINuE 
R2SCR  =  0.  DO 
DO  99  J  =1,3 

pjsGR)=~E2iQfiJi  P  ( I J  J  |  *  *  2 
CONTINUE 

ADJUST  ORIGINAL  RAY  TRANSIT  TIME 
1ST  (I)  =  LST  (I)  *  (DS  QRT  (R2S2R/R1SQR)  ) 

CONTINUE 

RETURN 


P  d,J) 


nnn  nnno  n  onnnno  n  nnnnnnnnnnnnnnn  no 


S3SSS3SSSSSSSSSSSSSSSSSSSSSSSSSS5SSS3SSSSSSSSSSSSSSSSSS3S 


SUBROUTINE  POSML?  (T 

THIS  FORTRAN  SUBROUTI 
IIKELIHOOD  PLANAR  MET 
APPARENT  POSITIONS. 

INPUTS  ARE  THE  HYDRO? 
POSITIONS,  AND  THE  YE 

OUTPUTS  ARE  THE  APPA5 
RAY  TRANSIT  TIMES  MT 
ALSO  OUTPUT  IS  A  VECT 
TIMING  ERROR  VARIANCE 
TO  THE  ACOUSTIC  CENTE 
COMPONENT  MEASURED  UP 

INTEGER  M,  I,  J 
DOUBLE  PRECISION  T  (13 
DOUBLE  PRECISION  V,TC 
DOUBLE  PRECISION  VAR( 

D  =  30. DO 
TOL  =  1.D-4 

INITIALIZE  THE  DIRECT 
USING  THE  LEAST  SQUAR 


,V,M,?,MT,  VAR) 

NE  IMPLEMENTS  THE  MAXIMUM 
HOD  M.L.P.  FOR  ESTIMATING 


HONE  TIMES  T  FOR  M  SOURCE 
LOCITY  OF  SOUND  V  AT  THE  ARRAY. 

ENT  POSITION  ESTIMATES  P  AND 
FOR  EACH  OF  THE  M  POSITIONS. 

OR  VAR  OF  M  ESTIMATES  OF  THE 
.  POSITIONS  ARE  ALL  1  EFEF.ENCED 
R  OF  THE  ARRAY,  WITH  THE  Z 
WARDS  FROM  THE  ARRAY. 


00,4)  ,?(  1000,3)  ,  MI  (10)  0) 
,D,DIFF,NUMER,jENOM,T)  L,C  (3) 
1000)  ,  CO  (3)  ,  X  (3) 


ION  COSINE  ESTIMATE  3Y 
ES  SOLUTION. 


LOOP  THROUGH  THE  M  POSITIONS 


DO  11  I  =  1 , M 
TC  =  T  ( I,  4) 

DENOM  =  (TC-T  (1,1) 

DO  12  J  =  1,3 


OM  =  (TC-T  (I,  1)  )  **2  +  {TC-T  (I,  2)  )  **2 

*  (TC-T  (1,3  **2 

12  J  =  1,3 

C(J)  =  (TC  -  I  ( I,  J) )  /D5QRT  (DENOM) 
TINUE 


CONTINUE 

ITER  =  0 
DENOM  =  0.  DO 
ITER  =  ITER  ♦  1 
TC  =  0. DO 

USE  ITERATION  FORM 
AND  TIME  VALUES. 

DO  33  J  =  1,3 
CO  (J)  =  C  (J) 
DENOM  =  DENOM  ♦ 
TC  =  TC  +  T  (I,J 
CONTINUE 

TC  =  (  (TC  +  T  (1,4)  )♦ 
DENOM  =  DENOM  -  TC 
DO  66  J  =  1,3 

xjjJ  =  iiAis*(cJ)  ( 

CONTINUE 


ULAE  TO  DEVELOPS  COSINES 


T(I,  J)*C(J) 


D*(C(1)  *CJ2)  ♦C(3))/V)  '4.D0 
*(C(T)*C(2)*C(j)J 


TC)  /DENOM 
J)-C(J)) 


CHECK  ITER  ATICN  TOLERANCES 

DIFF  =  DMAX1  (X  ( 1)  -X(2),X  (3)) 
IF  I  ITER  .LT.  100)  GO  TO  70 
FORMAT  (2X, 'ITERATIONS  EX 
FORMAT  ( 2X,  '  IN  POSITION 
WR I T  E  (6  , 60)  ITER 
WRITE  6,61)  I 
GD  TO  11 


POSITION  NO. 
R 


EXCEED  ',16) 


IF  (  TOL  -LI.  DIFF  )  GO  TO  22 

NU  MER  =  0.D0 
DENOM  =  0.  DO 


ACOUSTIC  C 


SINES  TO  POSITIONS,  AND  TRAMS.  AT] 
ENTER  REFERENCED  COORDINATES. 


DO  44  J  =  1,2 

P  (I,  J)  =  V* C  (J)  *TC 
DE  Nu  M  =  DENOM  *  P  (I,  J)  *  *2 
=  P(I,J)  -T5.D0 


ni(m£r  = 

CONTINUE 
NT  (I)  =  TC*DSQRT (NUMER/DENOM) 


NUMER  +  P  (I, J) **2 


ESTIMATE  VARIANCE 


VAR  (I)  =  0.  DO 
DO  55  J  =  1,3 

VAR  (I)  =  VAR  (I)  +  (T  (I,  J)  -TC+D*C(J)  /I)  **2 

CONTINUE 

VAR  (I )  =  VAR(I)/4.  DO 
CONTINUE 
RETURN 
END 


non  n  non  non  nnnnnnnn  no 


5  sssssssssss  ssss  ssssssssssssssssssssssssssssssssssssssssss 

SUEEOUTINE  POSMLS  (TIKE  ,V  ,M,  P,  MLT,  VARML) 


THIS  FORTRAN  SUBROUTINE  IMPLEMENTS  THE  MAXC  HUM 
LIKELIHOOD  SPHERICAL  METHOD  M.L.3.  FOR  ESTIMATING 
APPARENT  POSITIONS. 


INPUTS  AND  OUTPUTS  ARE  THE  SAME  AS  FOR  THE 
SUBROUTINE  POSMIP. 


SET  GIOEAL  VALUES 


R  =  (  3.  D  0- DSQR  T  (5 .  DO)  )/2.D0 
D  =  30.  DO 


TO  L 1  =  1.D-5 
TOI2  =  1.D-5 


START  OUTERMOST  LOOP  (ONE  FOR  EACH  SOURCE  PO>  ITIGN) 


DO  5  II  =  1,M 
DENCM  =  0. DO 
TAD  =  TI ME  (11,4) 
EPS  =  1.D-3 
TC  =  TAU 
DO  111  J  =  1,3 


11  1 


T(J)  =  TlfiE(II,J) 
Cjj  =  (  (  V**2)  *  { T1 
DENOM  =  DENOM  +  C 


122 


CONTINUE 

DO  122 
C 

CONT 

START  MIDDLE  LOOP 


C**2-T  (J)  **2)  /  (D*2.  D)  ))  ♦D/2 .DO 
(J)  **2 


cm  J=”c(5)3/dsqrt  (deno'm) 

TINUE 


(CREATE  DERIVATIVE  MATRIX  G?  (  ,)  ) 


10 


310 

31  1 


ITER  1  =  0 
II  =  TAU 
ITER  1  =  ITER  1 
IF  ]ITER1  . LT 


300 


OR  M A  T (2  X , 'EXC 
FORMAIf2X,’ 

TE  {6,3  101 

■■'I'3’1 


+  1 

50  )  GO  TO  300 

ESSIVE  NEWTON  ITERS. 
IN  POSITION 


WRI 
WRITE 
GO  TO 
L  =  0.  DO 
S  =  0. DO 
DLDT  =  0  .DO 
DSDT  =  0 . DO 


,16) 


II 


DO  11  I  =  1,3 

Kcff*  =~TAli*)*2*  (D/V) 
*<?)  =  WfltfDSQk 


VI 

IKill 


DLDC 

CSDC 


ill  : 

DLDT  =  DLD 


**2- (  2.  D0*D*TAU*C  ([  )  )  /V 


=  ({V*K 


m 


5) 


C(I)  *T  (I)  *D*  TAU)  /V  K  ( I) 

(I) 


C(I)  *T  (I)  *  (D*C  (I)  -V*TAU)  )  /VK  (I) 


nonn  nnn  non  nnn 


DSDT  =  DSDT+T  (I)*  D*C  (I)  -V*TAL 
CONTINUE 


-V  *TAU)  /VK  (I) 


DO  22  I  =  1,3 
DO  33  J  =  1,3 

GP  (I,J)  =  (TK(I)  *DLDC  (J)  )  /(L*L*DSQP.T(  K  (I)  )  ) 
CONTINUE 


CONTINUE 

GP  (4,1)  =  D SDC  (I)  *  (V*TC-D*L) 

GP  (4,1)  =  (GP  (4,1)  -D*DLDC  (I)  ♦  (1  -  DO— S)  ) 

/(7*(1.  DO-S)  *  *2 ) 

G?  (I,  I)  =  (1*T  (I)  *  D*TAU1  -DLDC(I)  *K  (I)  *V*TK  (I) 

GP  1,1)  =  1 .  DO  -  (GP  (1,1)  /  (VK  (I)  *L*L)  i 

GP  1,4)  =  (T  (I)  *L*  (V  *TAU -D*C  (I)  )  ) 

+  V  *K  (I)  *r  K  (I)  * 


GP  (1 , 4)  =  GP  (1,4)  /  (VK(I)*L*L) 

CONTINUE 

GP  (4,4)  =  1.D0-(  (DSDT  *(V*TC-D*L))-D*DLD?*M.  DO-S)  ) 

'  /  (V*  (1.  DO-S)  ^ ^2) 

INVEST  DERIVATIVE  MATRIX 

CALI  GAUSS3  (4,0,  GP,  GPI,  IER,  4) 

CALCULATE  INITIAL  NEWTON  SEARCH  SOLUTION 

DO  44  I  =1,3 

GO  (I)  =  C  (I)-TK  (I)/ (DSQRT(K  (I)  )  *L) 

CONTINUE 


*K  (I)  *C  K  (I)  *DLDT 


CONTIf 


4GC(4^Ajj  T  AD-  ( V*TC-D  *L)  /  ( V*  (1.D0-S)  ) 

D04^5  I  =U1,3 

lili  =  B  (4) -GPI  (4,1)  *GC(I) 

DO  56  J  =  1,4 


DO  56  J  =  1,4 

B  (I )  =  E  (I) -GPI  (I,J)  *GC(J) 
CONTINUE 
CONTINUE 

B  ( 4)  =  B  (4)  -  GPI  (4  ,4)  *GC  (4) 


PREPARE  FOR  GOLDEN  SECTION  SEARCH 


EPS  = 
DO  66 


=  DMAX1  ((1.D-6)  ,  (EPS/1. D1)) 
66  X  “  1,3 

*b*  *  b<i»  -* (i)  * 


CONTINUE 
A  (4)  =  TAU 

XT  (4)  =  A  (4)  +E*  (B  (4  )  -A  (4)  ) 

TAt)  =  XI  (4) 

START  INNERMOST  LOOP  {  COMMENCE  GOLDEN  SECTION  SEARCH 
TO  IMPROVE  THE  INITIAL  NEWTON  SOLUTION  ) 

ITER2  =  0 

CALL  OBJFCN  (FI  ,  L,  K  , TK  ,S  ,  XI  , D  ,  V,  T,  TC) 

ITER2  =  ITER 2  ♦  1 

IF  (ITER  2  .  LT.  50  )  GO  TO  4  00 

FORHAT(2X, 'EXCESSIVE  G.S.S.  ITERS.,  P)  S  #  *  ,16) 
WRITE  (6,4  10)  II 
GO  TO  5 
DIFF  =  0 .DO 
DO  77  I  =  1,3 

DIFF  =  DIFF  ♦  (A  (I)  -B  (I)  )  **2 

Stf'.'iilS  '  B'lj  '  X,(I) 

TINU] 


CC  NT 


nonn  noon  n  non  n  non 


33 

30 

99 

100 
144 


155 

5 

C 

csoo 


x2j4!  =  _A(4)  +  B  (4)  -  XI  (4) 


Ill 


X2 


DIFF  =  D^It  (DIFF  +  (A  (4) -B(4)  )  **2) 


CHECK  TOLERANCES  -  STOP  G-S.S.  IF  TIGHT  EXOJ  GH 
IF  (  EPS  .GT.  DIFF  )  GO  TO  100 
CALL  OBJFCN  (F2,L,K,  TK,  S, X2,D,  V,T,TC) 

CHOOSE  IH  PROVING  DIRECTION 


IF 

DO 


FI  .LT.  F2  )  30  TO  90 

81  =1.3 


CCNIxn!iE  X11 
A  (41  =  X  1(41 
XT  (4)  =  X2  (4) 
GO  TO  7  0 


DO  99  I  =  1 

w- 


=  x3  m 

)  =.  A  (IWB(I)-XI(I) 


: (I)  =  xi  ,  . 

CONTINUE 
B  ( 4)  =  X2  (4 ) 

XT  (4)  ir  A  (4  j  +B  (4)  -X  1  (4) 
GO  TO  70 


CHECK  TOLERANCES  FOR  NEWTON  SEARCH  ITERATIONS, 

AND  PREP  FOR  NEXT  NEWTON  ITERATION  IF  NECESSARY 


=  1,3 


DO  144  I  - 

X(I)  =  DArS  (C  (I)  -Cl  {I)) 
fTINUE 


CO  NT] 


DIFF  »  DMAX1  (X  (1)  ,X  (2)  ,X  (3)  ) 

IF  (  TO  L 1  .LT.  DIFF  )  GO  TO  10 
DIFF  =  DABS  {TAU-T 1 ^ 


DABS  (TAU 
IF  (  TOL2  .LI.  DIF 


)  GO  TO  10 


DONE  WITH  II-TH  SET  OF  TINES  AND  POSITIONS.  SAKE 
ESTIMATES,  AND  GC  ON  TO  NEXT  POSITION  TO  3E  ESTIMATED 


VAR  ML (I1 1  = 
DENOM  =  0. DO 


=  0.  DO 


NUMER  =  0. DO 
DO  155  J  =  1,3 


VAR  ML  (II )  =  VARHL(II)  ♦  If 
PJII.J)  =  V*TAU*C  J) 

DENOM  =  DENOM  ♦  P(II,J)**] 


IK  (J)  *TK  (J) 


(II,  J)  * *2 

P(II,J)  =  V*TAU*C  (J1-D/2.D0 
NUMER  =  NUMER  ♦  P(II,J)**2 
CONTINUE 


VARML(II)  =  (VARML  ( II)  +  (TC-TA  U)  * 
SDEV(II)  =  DSCRT  (VA  RMx.  (II)  1 
MLT(II)  =  TAE*DSQBT  (NUMER/DENOM) 
CONTINUE 

WRITER. 900)  (SDSV{I),I  =  1,M) 
FORMAT  ($E15. 5) 


*2) /4.D0 


RETURN 
END 


nnnnnnn  nn 


SSSSSSSSSSSSS  SSS SSSSSSS3S SSSS3SSSSSS3SS3SS5SSS>  SSSSSSSSSSS 

SUBROUTINE  OBJFCN  (F,I,K,TK,S,X,D,V,T,TC) 

THIS  FORTRAN  SUBROUTINE  IS  CALLED  BY  THE  SJBRCUTINE 
POSMLS.  IT  CALCULATES  NEW  VALUES  FOE  SEVERAL 
VARIABLES  AS  WELL  AS  A  FUNCTIONAL  VALUE  WHICH  i.S  ih’r. 
DECISION  FACTOF.  DETERMINING  THE  APPROPRIATE 
I M PROVING  DIRECTION  FOR  THE  GOLDEN  SECTION  SEARCH. 


DOUBLE  PRECISION  F,  L  ,  K(3)  ,  T  K  { 3)  ,  S , X ( 4)  ,DtV,  TC  ,  T  (3) 
C 

S  =  O.DO 
L  =  0.D0 
DO  11  I 

T^i) 

s  =  s 

L  =  L 
11  CONTINUE 
F  =  0  .  DO 
DO  22  I  =  1,3 

F  =  F+  (X(I)  -TK  (I)  /  (DSQRT  (K  (I)  )  *L)  )  **2 
22  CONTINUE 

F  =  F+  (X  (4)- (V*TC-D*L)/ {V*  (1.  DO-S)  )  )  **2 
RETURN 
END 


=  1,3 

=  X  (4 )  **2+  (D/VI  **2-  (  ( 2.  DQ*D*  X  (4)  *X  ( 
=  T (I) -DSQRT  (K  (I) ) 

♦  TK  (I)  /DSQRT  (K  (lj  ) 

+  (X  I)  *TX(I)  )/D~  '  '  * 


:))/v) 


*TX 


SQivI  (K  (I)  ) 


114 
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